---
title: "Get started with EpiExposure"
author: "Ricardo G. Tomáz"
description: "A first exposure-lag analysis with EpiExposure"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Get started with EpiExposure}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

library(EpiExposure)
library(ggplot2)
library(dplyr)
```

# What does `EpiExposure` do?

`EpiExposure` provides a standardized workflow for constructing, fitting, interpreting, and simulating distributed exposure-lag models in plant disease epidemiology.

The package helps users:

- define non-linear exposure-lag structures;
- fit distributed lag non-linear models;
- quantify lag-specific, period-specific, and cumulative effects;
- evaluate complete exposure profiles;
- predict expected disease outcomes;
- simulate alternative exposure scenarios;
- estimate yield and economic losses.

Unlike conventional window-based approaches, `EpiExposure` treats the complete environmental profile preceding disease assessment as the epidemiological exposure unit.

::: tutorial-note
**What you will learn.** This tutorial shows how to organize
longitudinal epidemiological data, define and diagnose an exposure–lag
structure, fit a first DLNM, and predict an expected disease outcome
from a complete exposure profile.
:::

# The EpiExposure workflow

A basic analysis follows six stages:

1. organize and inspect complete exposure profiles;
2. define the exposure-response and lag-response basis functions;
3. evaluate the proposed structure for identifiability and numerical
   stability;
4. collapse each profile into one epidemic-level design row;
5. prepare the response and fit the model;
6. use the fitted model for interpretation, prediction, or simulation.

This section introduces the basic workflow but stops before detailed effect interpretation. Lag-specific, period-specific, and cumulative effects are covered in the *Understanding Effects* section

# Organizing exposure profiles

`EpiExposure` expects longitudinal data in which each row represents one observation within an epidemic profile.

At minimum, the dataset must contain:

1.  an epidemic identifier, such as `epi_id`;
2.  a numeric temporal index, such as `time`;
3.  one or more exposure variables;
4.  one final response associated with each epidemic.

In the example dataset:

- `epi_id` identifies independent epidemic profiles;
- `time` orders observations chronologically;
- `tmean`, `rain`, and `wetness` are environmental exposures;
- `y` is the final disease response.

The temporal index begins at the earliest observation and ends at the disease assessment:

``` text
time 0  → earliest exposure observation → maximum lag
time 85 → disease assessment            → lag 0
```

Every epidemic must contain the same complete and regularly spaced exposure profile. With `max_lag = 85`, each epidemic requires:

$$85 + 1 = 86$$

observations for every fitted exposure.

```{r load-package}
library(EpiExposure)
library(ggplot2)
```

```{r load-data}
data("epi_data")
```

# Inspecting the example dataset

The example contains complete environmental profiles for multiple epidemics.

```{r example-data-table, echo=FALSE}
knitr::kable(
  utils::head(
    epi_data,
    10L
  ),
  digits = 2,
  align = "c"
)
```

Before modeling, verify that every epidemic contains the required number of observations:

```{r check-history-length}
table(
  table(
    epi_data$epi_id
  )
)
```

For `max_lag = 85`, the expected result is 86 observations per `epi_id`.

# Visualizing exposure profiles

Visual inspection can reveal temporal trends, unusual exposure profiles, limited variation, and missing observations.

```{r, message=FALSE, warning=FALSE}
wetness_plot = epi_data |>  
  ggplot(aes(time,wetness, group = epi_id))+
  geom_smooth(se = F, color = "#55A868",linewidth = 0.3)+
  theme_bw()+
  labs(x = "Days after planting the crop",
       y = "Leaf wetness (%)")+
  scale_x_continuous(limits = c(0,85), expand = c(0,0),
                     breaks = pretty(epi_data$time, n = 8))+
  theme(text = element_text(size = 10))


wetness_plot
```


```{r, message=FALSE, warning=FALSE}
rain_plot = epi_data  |>  
  ggplot(aes(time,rain, group = epi_id))+
  geom_smooth(se = F, color = "#4C72B0",linewidth = 0.3)+
  theme_bw()+
  labs(x = "Days after planting the crop",
       y = "Daily precipitation (mm)")+
  scale_x_continuous(limits = c(0,85), expand = c(0,0),
                     breaks = pretty(epi_data$time, n = 8))+
  theme(text = element_text(size = 10))


rain_plot
```


```{r, message=FALSE, warning=FALSE}
tmean_plot = epi_data |> 
  ggplot(aes(time,tmean, group = epi_id))+
  geom_smooth(se = F, color = "#C44E52",linewidth = 0.3)+
  theme_bw()+
  labs(x = "Days after planting the crop",
       y = "Mean temperature (°C)")+
  scale_x_continuous(limits = c(0,85), expand = c(0,0),
                     breaks = pretty(epi_data$time, n = 8))+
  theme(text = element_text(size = 10))


tmean_plot
```
The objective is not to identify exposure effects visually. These plots only describe the observed profiles. Effect estimation requires the DLNM model.

# Defining exposure-lag structures

The first modeling step is to define the non-linear exposure-response and lag-response basis functions.

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

Here:

- `vars` identifies the exposures included in the DLNM;
- `max_lag = 85` defines the retrospective lag window;
- `df_var = 3` controls exposure-response flexibility;
- `df_lag = 2` controls lag-response flexibility.

The resulting object stores the fitted basis definitions required for model construction and downstream prediction.

::: tutorial-note
**Exposure variables only.** Coordinates, treatment groups, blocks, years, and other identifiers should not be included in `vars`. They can be preserved separately by `build_design()` when needed for random or spatial effects.
:::

# Checking identifiability and numerical stability

Before fitting the model, the proposed exposure–lag structure should be
evaluated for identifiability and numerical stability. Increasing the
number of exposures or the flexibility of the exposure and lag bases
increases the number of cross-basis predictors and may produce
redundant, rank-deficient, or poorly conditioned designs.

`check_identifiability()` reconstructs the same epidemic-level
cross-basis design used by the EpiExposure fitting workflow. The
function evaluates the combined numerical rank, scaled condition number,
near-zero-variance columns, between-exposure cross-basis correlations,
supplementary variance inflation factors, and the number of epidemics
relative to design complexity.

```{r check-identifiability}
identifiability <- check_identifiability(
  data = epi_data,
  vars = c(
    "tmean",
    "rain",
    "wetness"
  ),
  max_lag = 85,
  df_var = 3,
  df_lag = 2
)


```

The printed output classifies the proposed design as `"ok"`, `"warning"`,
or `"problem"`. The `identifiable` field indicates whether the combined
epidemic-level design has full numerical rank, whereas
`numerically_stable` additionally considers severe conditioning,
near-zero-variance columns, and whether the available epidemics can
support the number of design columns.

```{r inspect-identifiability1}
knitr::kable(
  utils::head(
    identifiability$overall,
    10L
  ),
  digits = 4,
  align = "c"
)
```

```{r inspect-identifiability2}
knitr::kable(
  utils::head(
    identifiability$by_variable,
    10L
  ),
  digits = 4,
  align = "c"
)
```

```{r inspect-identifiability3}
knitr::kable(
  utils::head(
    identifiability$recommendations,
    10L
  ),
  digits = 4,
  align = "c"
)
```

A `"warning"` does not automatically invalidate the proposed structure.
It indicates that the reported diagnostics should be examined before
model fitting. A `"problem"` indicates rank deficiency or severe
numerical instability and generally supports simplifying the candidate
structure, for example by reducing `df_var` or `df_lag`, reconsidering
redundant exposures, or increasing the number of independent epidemics.

The diagnostic thresholds are interpretive heuristics rather than
universal inferential rules. In particular, high correlations or VIFs
among individual spline columns may arise from the basis construction
itself. The combined numerical rank and scaled condition number should
therefore receive greater emphasis than any isolated column-level
diagnostic.

# Building the epidemic-level design

`build_design()` transforms each complete longitudinal profile into one epidemic-level row containing the cross-basis predictors.

```{r build-design}
dat <- build_design(
  data = epi_data,
  cb_templates = cb)
```

The original exposure columns are represented by canonical cross-basis columns such as:

``` text
cb_tmean_1, cb_tmean_2, ...
cb_rain_1, cb_rain_2, ...
cb_wetness_1, cb_wetness_2, ...
```

The coefficients of these columns should not be interpreted individually. Together, they reconstruct the non-linear exposure-lag-response association.

Inspect the resulting design:

```{r design-table, echo=FALSE}
knitr::kable(
  utils::head(
    dat,
    6L
  ),
  digits = 3,
  caption = "First six rows of the epidemic-level DLNM design.",
  align = "c"
)
```

# Preparing the disease response

The response must be prepared for the distribution used in the model.

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

The Beta family is appropriate when the response is a continuous proportion strictly between 0 and 1.

Other supported families should be selected according to the response scale, for example:

- `binomial` for a Bernoulli response coded as 0 or 1;
- `poisson` for non-negative counts;
- `negative_binomial` for overdispersed counts;
- `gamma` for positive continuous responses;
- `gaussian` for approximately continuous, unbounded responses.

For details on available modeling engines, supported response distributions, engine-specific requirements, and Bayesian and frequentist workflows, see the *Advanced Topics* section.

# Fitting the first DLNM

The prepared design can now be fitted with one of the supported modeling engines.

```{r fit-model}
fit <- fit_epidlnm(
  data = dat,
  model_engine = "glmmTMB",
  family = "beta",
  random_effect = NULL,
  epiexposure_spec = attr(
    cb,
    "spec"
  )
)
```

This model estimates a population-level association between the complete profiles of temperature, rainfall, and leaf wetness and the final disease response.

```{r model-summary}
summary(fit)
```

The individual `cb_*` coefficients are components of the cross-basis expansion. Their signs and $p$-values should not be interpreted as isolated epidemiological effects.

The fitted exposure-lag surface must instead be reconstructed into quantities such as lag-specific, period-specific, and cumulative effects. That process is introduced in the *Understanding Effects* section

# Predicting one complete exposure profile

A fitted model can also estimate the expected disease outcome associated with a complete hypothetical exposure profile.

In this example, a constant environmental scenario is evaluated over an 86-day exposure window. Daily mean temperature (`tmean`) is fixed at 25°C, daily rainfall (`rain`) at 5 mm, and daily leaf wetness duration (`wetness`) at 15 hours. The model combines the lagged contributions of all daily exposures to estimate the expected disease outcome for this exposure profile.

```{r first-prediction}
prediction <- predict_outcomes(
  fit = fit,
  profiles = list(
    tmean = rep(
      25,
      86
    ),
    rain = rep(
      5,
      86
    ),
    wetness = rep(
      15,
      86
    )
  )
)

```

```{r}
prediction
```

Each exposure profile contains exactly 86 chronologically ordered observations, matching:

``` text
max_lag + 1
```

The resulting prediction represents the expected population-level disease response for the complete exposure profile. It is not a newly simulated observation.

# Where to go next

You have now completed the basic `EpiExposure` workflow:

1. loaded longitudinal epidemic data;
2. inspected complete exposure profiles;
3. defined the exposure–lag basis;
4. evaluated its identifiability and numerical stability;
5. constructed the epidemic-level design;
6. prepared the response;
7. fitted a DLNM;
8. predicted an expected disease outcome.

Continue with the *Understanding Effects* section to learn how to:

- interpret lag-specific effects;
- define epidemiological periods;
- calculate cumulative effects;
- select meaningful reference exposures;
- visualize uncertainty across the exposure-lag surface.

More advanced topics, including alternative modeling engines, random
effects, spatial Matérn structures, and Bayesian estimation, are covered
in the *Advanced Topics* section. Additional comparisons among candidate
exposure–lag structures are discussed in the *Sensitivity and Decision*
section.