Introduction

Understanding exposure-lag associations is important for epidemiological interpretation, but researchers are often interested in a different question:

What is the expected disease outcome under an alternative environmental profile?

EpiExposure provides simulation and prediction tools that allow researchers to:

• generate complete synthetic exposure profiles;
• introduce structured exposure events;
• construct hypothetical epidemiological scenarios;
• predict expected disease outcomes;
• propagate parameter uncertainty;
• compare alternative exposure profiles and their predictions.

Unlike lag-specific effect summaries, prediction evaluates all exposure values and lags jointly. Therefore, every prediction represents a complete environmental profile rather than an isolated exposure value.

The outputs describe predictions and contrasts under the fitted statistical model. Hypothetical scenarios should not automatically be interpreted as causal interventions unless the study design and assumptions support that interpretation.

Preparing the example model

This section uses the same dataset and DLNM specification introduced in the previous tutorials.

library(EpiExposure)
library(ggplot2)

data("epi_data")

Define the exposure-lag structures:

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

Build the epidemic-level design:

dat <- build_design(
  data = epi_data,
  cb_templates = cb
)

Prepare the response:

dat <- prepare_response(
  data = dat,
  response = "y",
  family = "beta"
)

Fit the model:

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

The fitted model relates complete profiles of temperature, rainfall, and leaf wetness to the expected disease response.

Complete exposure profiles

Each exposure profile supplied to predict_outcomes() must contain exactly max_lag + 1 observations. With max_lag = 85, each profile therefore contains 86 values.

Profiles are supplied in chronological order, from the earliest exposure observation to the most recent:

• the first value (time = 0) is the oldest observation and corresponds to lag 85;
• the final value (time = 85) is the most recent observation and corresponds to lag 0.

More generally, chronological time and retrospective lag are related by:

\[ \text{lag} = \text{max\_lag} - \text{time}. \]

Profile construction and visualization use chronological time, whereas DLNM effects are interpreted on the retrospective lag scale. Thus, chronological time increases from the oldest to the most recent observation, while retrospective lag increases in the opposite direction, from the most recent to progressively older observations.

simulate_exposures() generates profiles directly in chronological order and does not reverse them. Simulated or manually defined profiles should therefore be supplied to predict_outcomes() in the same order in which they were created. Users should not reverse the profiles manually.

Simulating exposure profiles

Simulating a background process

The following example generates a temperature profile using an autoregressive background process:

temperature_profile <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 123,
  background = list(
    dist = "ar1",
    mean = 25,
    sd = 4,
    phi = 0.90
  )
)

The result contains the simulated chronological profile:

#>   time_0   time_1   time_2   time_3   time_4   time_5   time_6   time_7 
#> 22.75810 22.58096 25.54056 25.60944 25.77392 28.68685 29.12180 26.50391 
#>   time_8   time_9  time_10  time_11  time_12  time_13  time_14  time_15 
#> 25.15595 24.36331 26.56124 27.03248 27.52800 27.46818 26.25222 29.24259 
#>  time_16  time_17  time_18  time_19  time_20  time_21  time_22  time_23 
#> 29.68636 25.78881 26.93278 25.91517 23.96184 23.68560 22.02814 21.05446 
#>  time_24  time_25  time_26  time_27  time_28  time_29  time_30  time_31 
#> 20.35922 17.88245 20.05494 20.81686 19.25076 22.01179 23.05417 22.73428 
#>  time_32  time_33  time_34  time_35  time_36  time_37  time_38  time_39 
#> 24.52156 26.10048 27.42291 28.38130 29.00896 28.50012 27.61664 26.69160 
#>  time_40  time_41  time_42  time_43  time_44  time_45  time_46  time_47 
#> 25.31118 24.91755 22.71950 26.72925 28.66248 26.33803 25.50177 24.63795 
#>  time_48  time_49  time_50  time_51  time_52  time_53  time_54  time_55 
#> 26.03407 25.78531 26.14845 25.98383 25.81070 28.11587 27.41064 29.81363 
#>  time_56  time_57  time_58  time_59  time_60  time_61  time_62  time_63 
#> 26.63193 27.48804 27.45519 27.58617 27.98948 26.81470 26.05226 24.17109 
#>  time_64  time_65  time_66  time_67  time_68  time_69  time_70  time_71 
#> 22.38525 23.17595 24.13983 24.31826 25.99447 29.46946 28.16638 23.82356 
#>  time_72  time_73  time_74  time_75  time_76  time_77  time_78  time_79 
#> 25.69477 24.38876 23.25030 25.21342 24.69556 22.59761 23.15396 23.09640 
#>  time_80  time_81  time_82  time_83  time_84  time_85 
#> 23.29681 24.13889 23.57873 24.84437 24.47550 25.10643

The profile is ordered from the oldest to the most recent exposure observation. For this 86-value profile, the first value corresponds to lag 85 and the final value corresponds to lag 0. The profile remains in chronological order when it is later supplied to predict_outcomes().

Visualize the complete profile:

df_tmean <- data.frame(
  time = 0:85,
  value = as.numeric(temperature_profile$profile)
)

 df_tmean%>% 
  ggplot(aes(time,value))+
  geom_line(color = "#C44E52", size = 1.4)+
  theme_bw()+
  labs(x = "Days after planting the crop",
       y = "Exposure")+
  scale_x_continuous(
  limits = c(0, 85),
  breaks = seq(0, 85, by = 10)
)+
  theme(text = element_text(size = 10, face = "bold"))

The autoregressive parameter:

phi = 0.90

produces temporal dependence among consecutive observations. It does not represent an exposure-lag effect. The lag-response relationship remains determined by the fitted DLNM.

Simulating structured exposure events

A structured event can be superimposed on a background profile.

The following example constructs a rainfall pattern with values assigned to selected portions of the temporal profile:

rain_profile <- simulate_exposures(
  max_lag = 85,
  mode = "pattern",
  pattern = "block_times",
  block_times = c(
    0,
    40
  ),
  fixed_value = 15,
  background = list(
    dist = "fixed",
    value = 0
  )
)

#>  time_0  time_1  time_2  time_3  time_4  time_5  time_6  time_7  time_8  time_9 
#>      15      15      15      15      15      15      15      15      15      15 
#> time_10 time_11 time_12 time_13 time_14 time_15 time_16 time_17 time_18 time_19 
#>      15      15      15      15      15      15      15      15      15      15 
#> time_20 time_21 time_22 time_23 time_24 time_25 time_26 time_27 time_28 time_29 
#>      15      15      15      15      15      15      15      15      15      15 
#> time_30 time_31 time_32 time_33 time_34 time_35 time_36 time_37 time_38 time_39 
#>      15      15      15      15      15      15      15      15      15      15 
#> time_40 time_41 time_42 time_43 time_44 time_45 time_46 time_47 time_48 time_49 
#>      15       0       0       0       0       0       0       0       0       0 
#> time_50 time_51 time_52 time_53 time_54 time_55 time_56 time_57 time_58 time_59 
#>       0       0       0       0       0       0       0       0       0       0 
#> time_60 time_61 time_62 time_63 time_64 time_65 time_66 time_67 time_68 time_69 
#>       0       0       0       0       0       0       0       0       0       0 
#> time_70 time_71 time_72 time_73 time_74 time_75 time_76 time_77 time_78 time_79 
#>       0       0       0       0       0       0       0       0       0       0 
#> time_80 time_81 time_82 time_83 time_84 time_85 
#>       0       0       0       0       0       0

Visualize the rainfall profile:

df_rain <- data.frame(
  time = 0:85,
  value = as.numeric(rain_profile$profile)
)

 df_rain%>% 
  ggplot(aes(time,value))+
  geom_line(color = "#4C72B0", size = 1.4)+
  theme_bw()+
  labs(x = "Days after planting the crop",
       y = "Exposure")+
  scale_x_continuous(
  limits = c(0, 85),
  breaks = seq(0, 85, by = 10)
)+
  theme(text = element_text(size = 10, face = "bold"))

A leaf wetness profile can be generated independently:

wetness_profile <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 456,
  background = list(
    dist = "ar1",
    mean = 10,
    sd = 4,
    phi = 0.90
  )
)

#>    time_0    time_1    time_2    time_3    time_4    time_5    time_6    time_7 
#>  4.625914  6.247426  8.019056  5.795534  4.970456  4.908391  6.621729  7.396401 
#>    time_8    time_9   time_10   time_11   time_12   time_13   time_14   time_15 
#>  9.413140 10.471295  8.827395 11.230632 12.831472 15.432048 12.376714 15.534374 
#>   time_16   time_17   time_18   time_19   time_20   time_21   time_22   time_23 
#> 18.009388 17.884050 21.071020 22.645309 20.553278 16.503720 13.365585 13.392098 
#>   time_24   time_25   time_26   time_27   time_28   time_29   time_30   time_31 
#> 12.990405 14.669058 13.395137 12.483066 14.823142 12.441433 11.275305 10.112462 
#>   time_32   time_33   time_34   time_35   time_36   time_37   time_38   time_39 
#>  6.615987  7.470749  8.021169 11.384750 10.094474  9.840489  9.117208  9.137949 
#>   time_40   time_41   time_42   time_43   time_44   time_45   time_46   time_47 
#>  9.173692  9.941607  9.512229  9.706507 13.360498 13.235160 13.117645 14.148516 
#>   time_48   time_49   time_50   time_51   time_52   time_53   time_54   time_55 
#> 11.684280 12.229035 10.846750 10.314809 11.466825 12.883845 13.673602 14.581575 
#>   time_56   time_57   time_58   time_59   time_60   time_61   time_62   time_63 
#> 13.403023 15.778860 16.145459 16.537209 11.952268 10.044761  9.690582 11.283632 
#>   time_64   time_65   time_66   time_67   time_68   time_69   time_70   time_71 
#> 11.079345 13.843590 16.184937 13.070750 11.814142 10.288002  9.962817  9.895465 
#>   time_72   time_73   time_74   time_75   time_76   time_77   time_78   time_79 
#>  6.906054  9.733334  9.725104 10.223759 10.004476  8.513481  6.859744  7.560886 
#>   time_80   time_81   time_82   time_83   time_84   time_85 
#>  5.922656  5.560829  9.527843  9.702088 10.967745 12.353884

df_wetness <- data.frame(
  time = 0:85,
  value = as.numeric(wetness_profile$profile)
)


 df_wetness%>% 
  ggplot(aes(time,value))+
  geom_line(color = "#55A868", size = 1.4)+
  theme_bw()+
  labs(x = "Days after planting the crop",
       y = "Exposure")+
  scale_x_continuous(
  limits = c(0, 85),
  breaks = seq(0, 85, by = 10)
)+
  theme(text = element_text(size = 10, face = "bold"))

The three profiles can now be combined into one complete hypothetical environmental profile.

All three profiles follow the same chronological convention and contain 86 observations. At each chronological position, the temperature, rainfall, and leaf wetness values describe the environmental conditions assigned to that point in the hypothetical profile.

When these profiles are supplied to predict_outcomes(), the function validates their length and order and projects each complete profile through the cross-basis definition stored in the fitted model. The resulting design follows the retrospective DLNM lag convention, in which lag 0 represents the most recent exposure and increasing lags represent progressively older exposures.

This internal transformation does not require the user to reverse or otherwise rearrange the simulated profiles.

Predicting disease outcomes

prediction <- predict_outcomes(
  fit = fit,
  profiles = list(
    tmean = temperature_profile$profile,
    rain = rain_profile$profile,
    wetness = wetness_profile$profile
  ),
  uncertainty = FALSE
)

prediction
#>   prediction
#> 1  0.1195241

Prediction uncertainty

A deterministic prediction uses the central fitted coefficient estimates. However, the coefficients are estimated with uncertainty.

When:

uncertainty = TRUE

EpiExposure propagates joint coefficient uncertainty through the complete prediction calculation.

To inspect individual uncertainty realizations, request sample-level output:

prediction_samples <- predict_outcomes(
  fit = fit,
  profiles = list(
    tmean = temperature_profile$profile,
    rain = rain_profile$profile,
    wetness = wetness_profile$profile
  ),
  uncertainty = TRUE,
  n_samples = 100,
  output = "samples",
  seed = 123
)

head(prediction_samples)
#>   sample prediction
#> 1      1  0.1034903
#> 2      2  0.1042327
#> 3      3  0.0892574
#> 4      4  0.1337358
#> 5      5  0.0837315
#> 6      6  0.1593514

Inspect the first realizations:

Visualize the prediction distribution:

ggplot(
  prediction_samples,
  aes(
    x = prediction * 100
  )
) +
  geom_histogram(
    bins = 25,
    fill = "#4C72B0",
    color = "white"
  ) +
  theme_bw() +
  labs(
    x = "Predicted disease response (%)",
    y = "Number of uncertainty realizations",
    title = "Parameter uncertainty in the predicted response"
  )

This distribution represents parameter uncertainty in the expected response. It does not include observation-level noise and should not be interpreted as a distribution of newly simulated disease observations.

A summarized prediction can also be requested when individual realizations are not needed:

prediction_summary <- predict_outcomes(
  fit = fit,
  profiles = list(
    tmean = temperature_profile$profile,
    rain = rain_profile$profile,
    wetness = wetness_profile$profile
  ),
  uncertainty = TRUE,
  n_samples = 100,
  output = "summary",
  seed = 123
)

prediction_summary
#>   prediction prediction_sd prediction_lower prediction_upper
#> 1  0.1293006    0.03879026       0.07271773        0.2212469

Scenario simulation

Defining epidemiological periods

The scenario grid can be organized using epidemiological periods:

periods <- define_periods(
  max_lag = 85,
  cuts = c(
    20,
    40,
    60
  ),
  prefix = "Period"
)

periods
#>    period lag_start lag_end
#> 1 Period1         0      20
#> 2 Period2        21      40
#> 3 Period3        41      60
#> 4 Period4        61      85

The resulting periods divide the retrospective lag interval into four subintervals.

Remember that lag order is retrospective:

smaller lag → closer to disease assessment
larger lag  → earlier exposure condition

Generating candidate scenarios

The following example varies rainfall across a grid while evaluating selected temperature and leaf wetness conditions:

scenarios <- simulate_ranges(
  periods = periods,
  vary = list(
    tmean = c(
      20,
      34
    ),
    rain = seq(
      0,
      15,
      by = 1
    ),
    wetness = c(
      2,
      10
    )
  ),
  scenario_type = "grid",
  scenario_var = "rain"
)

Predicting scenario outcomes

simulate_scenarios() evaluates the expected disease outcome associated with each complete exposure scenario. When uncertainty = TRUE, uncertainty in the joint fitted model parameters is propagated through the prediction calculation.

output = "summary"

pred_scenarios <- simulate_scenarios(
  fit = fit,
  scenarios = scenarios,
  data = epi_data,
  uncertainty = TRUE,
  output = "summary",
  n_samples = 100,
  seed = 123
)

plot_scenarios(
  data = pred_scenarios,
  x = "rain",
  facet = c(
    "wetness",
    "tmean"
  ),
  x_limits = c(
    0,
    15
  ),
  output = "summary"
)

output = "samples"

pred_scenarios_s <- simulate_scenarios(
  fit = fit,
  scenarios = scenarios,
  data = epi_data,
  uncertainty = TRUE,
  output = "samples",
  n_samples = 100,
  seed = 123
)

plot_scenarios(
  data = pred_scenarios_s,
  x = "rain",
  facet = c(
    "wetness",
    "tmean"
  ),
  x_limits = c(
    0,
    15
  ),
  output = "samples"
)

Scenario surfaces

When enough values are evaluated along two continuous dimensions, predictions can be presented as a scenario surface.

The current example evaluates only two temperature values. Therefore, a line or point representation is more appropriate than a dense continuous surface.

A tile display can still summarize the evaluated grid:

ggplot(
  pred_scenarios,
  aes(
    x = rain,
    y = factor(
      tmean
    ),
    fill = prediction * 100
  )
) +
  geom_tile() +
  facet_wrap(
    vars(
      wetness
    ),
    labeller = label_both
  ) +
  scale_fill_viridis_c(
    name = "Expected\ndisease (%)"
  ) +
  theme_bw() +
  labs(
    x = "Rainfall",
    y = "Mean temperature",
    title = "Predicted response across the evaluated scenario grid"
  )

A smoother surface requires a denser sequence of temperature values in simulate_ranges().

Comparing predictions

exposure_tmean_A <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 101,
  background = list(
    dist = "ar1",
    mean = 25,
    sd = 2,
    phi = 0.90
  ),
  n = 1
)

exposure_rain_A <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 102,
  background = list(
    dist = "ar1",
    mean = 5,
    sd = 2,
    phi = 0.50
  ),
  n = 1
)

exposure_wetness_A <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 103,
  background = list(
    dist = "ar1",
    mean = 15,
    sd = 3,
    phi = 0.70
  ),
  n = 1
)

scenario_A <- list(
  tmean = exposure_tmean_A,
  rain = exposure_rain_A,
  wetness = exposure_wetness_A
)

exposure_tmean_B <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 201,
  background = list(
    dist = "ar1",
    mean = 34,
    sd = 2,
    phi = 0.95
  ),
  n = 1
)

exposure_rain_B <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 202,
  background = list(
    dist = "ar1",
    mean = 5,
    sd = 1.5,
    phi = 0.40
  ),
  n = 1
)

exposure_wetness_B <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 203,
  background = list(
    dist = "ar1",
    mean = 2,
    sd = 0.7,
    phi = 0.60
  ),
  n = 1
)

scenario_B <- list(
  tmean = exposure_tmean_B,
  rain = exposure_rain_B,
  wetness = exposure_wetness_B
)

exposure_tmean_C <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 301,
  background = list(
    dist = "ar1",
    mean = 15,
    sd = 2,
    phi = 0.90
  ),
  n = 1
)

exposure_rain_C <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 302,
  background = list(
    dist = "ar1",
    mean = 0.5,
    sd = 0.5,
    phi = 0.40
  ),
  n = 1
)

exposure_wetness_C <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 303,
  background = list(
    dist = "ar1",
    mean = 2,
    sd = 0.8,
    phi = 0.60
  ),
  n = 1
)

scenario_C <- list(
  tmean = exposure_tmean_C,
  rain = exposure_rain_C,
  wetness = exposure_wetness_C
)

exposure_tmean_D <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 401,
  background = list(
    dist = "ar1",
    mean = 10,
    sd = 2,
    phi = 0.90
  ),
  n = 1
)

exposure_rain_D <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 402,
  background = list(
    dist = "ar1",
    mean = 15,
    sd = 4,
    phi = 0.60
  ),
  n = 1
)

exposure_wetness_D <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 403,
  background = list(
    dist = "ar1",
    mean = 20,
    sd = 3,
    phi = 0.70
  ),
  n = 1
)

scenario_D <- list(
  tmean = exposure_tmean_D,
  rain = exposure_rain_D,
  wetness = exposure_wetness_D
)

profile_to_data <- function(
    exposure,
    variable,
    scenario
) {
  data.frame(
    scenario = scenario,
    var = variable,
    time = 0:85,
    lag = 85:0,
    value = as.numeric(exposure$profile),
    stringsAsFactors = FALSE
  )
}

df_simulations <- dplyr::bind_rows(
  profile_to_data(
    exposure_tmean_A,
    variable = "tmean",
    scenario = "A"
  ),
  profile_to_data(
    exposure_rain_A,
    variable = "rain",
    scenario = "A"
  ),
  profile_to_data(
    exposure_wetness_A,
    variable = "wetness",
    scenario = "A"
  ),
  profile_to_data(
    exposure_tmean_B,
    variable = "tmean",
    scenario = "B"
  ),
  profile_to_data(
    exposure_rain_B,
    variable = "rain",
    scenario = "B"
  ),
  profile_to_data(
    exposure_wetness_B,
    variable = "wetness",
    scenario = "B"
  ),
  profile_to_data(
    exposure_tmean_C,
    variable = "tmean",
    scenario = "C"
  ),
  profile_to_data(
    exposure_rain_C,
    variable = "rain",
    scenario = "C"
  ),
  profile_to_data(
    exposure_wetness_C,
    variable = "wetness",
    scenario = "C"
  ),
  profile_to_data(
    exposure_tmean_D,
    variable = "tmean",
    scenario = "D"
  ),
  profile_to_data(
    exposure_rain_D,
    variable = "rain",
    scenario = "D"
  ),
  profile_to_data(
    exposure_wetness_D,
    variable = "wetness",
    scenario = "D"
  )
)

simulation_plot <- ggplot(
  df_simulations,
  aes(
    x = time,
    y = value,
    color = var
  )
) +
  geom_line(linewidth = 1) +
  facet_grid(
    rows = vars(var),
    cols = vars(scenario),
    scales = "free_y"
  ) +
  scale_x_continuous(
    limits = c(
      0,
      85
    ),
    breaks = seq(
      0,
      85,
      by = 15
    )
  ) +
  scale_color_manual(
    values = c(
      tmean = "#C44E52",
      rain = "#4C72B0",
      wetness = "#55A868"
    )
  ) +
  theme_bw() +
  labs(
    x = "Chronological time",
    y = "Exposure",
    color = "Variable"
  ) +
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    ),
    legend.position = "none",
    strip.background = element_rect(
      color = "black",
      fill = "white"
    ),
    strip.text = element_text(
      face = "bold"
    )
  )

simulation_plot

comparison_tmean_AD <- compare_exposures(
  exposure1 = exposure_tmean_A$profile,
  exposure2 = exposure_tmean_D$profile,
  mode = "timewise"
)

comparison_tmean_AD$var <- "tmean"
comparison_tmean_AD$scenario1 <- "A"
comparison_tmean_AD$scenario2 <- "D"

comparison_rain_AD <- compare_exposures(
  exposure1 = exposure_rain_A$profile,
  exposure2 = exposure_rain_D$profile,
  mode = "timewise"
)

comparison_rain_AD$var <- "rain"
comparison_rain_AD$scenario1 <- "A"
comparison_rain_AD$scenario2 <- "D"

comparison_wetness_AD <- compare_exposures(
  exposure1 = exposure_wetness_A$profile,
  exposure2 = exposure_wetness_D$profile,
  mode = "timewise"
)

comparison_wetness_AD$var <- "wetness"
comparison_wetness_AD$scenario1 <- "A"
comparison_wetness_AD$scenario2 <- "D"

exposure_comparison_AD <- dplyr::bind_rows(
  comparison_tmean_AD,
  comparison_rain_AD,
  comparison_wetness_AD
)

exposure_comparison_AD$var <- factor(
  exposure_comparison_AD$var,
  levels = c(
    "tmean",
    "rain",
    "wetness"
  )
)

exposure_comparison_plot <- ggplot(
  exposure_comparison_AD,
  aes(
    x = time,
    y = diff,
    color = var
  )
) +
  geom_hline(
    yintercept = 0,
    color = "grey50",
    linetype = "dashed",
    linewidth = 0.5
  ) +
  geom_line(
    linewidth = 1
  ) +
  facet_wrap(
    vars(var),
    ncol = 1,
    scales = "free_y",
    labeller = as_labeller(
      c(
        tmean = "Mean temperature",
        rain = "Rainfall",
        wetness = "Leaf wetness"
      )
    )
  ) +
  scale_x_continuous(
    limits = c(
      0,
      85
    ),
    breaks = seq(
      0,
      85,
      by = 15
    )
  ) +
  scale_color_manual(
    values = c(
      tmean = "#C44E52",
      rain = "#4C72B0",
      wetness = "#55A868"
    )
  ) +
  theme_bw() +
  labs(
    x = "Days after planting the crop",
    y = "Scenario D - Scenario A",
    color = "Exposure"
  ) +
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    ),
    legend.position = "none",
    strip.background = element_rect(
      fill = "white",
      color = "black"
    ),
    strip.text = element_text(
      face = "bold"
    )
  )

exposure_comparison_plot

output = "summary"

comp_pred <- compare_predictions(
  fit = fit,
  profiles = list(
    A = scenario_A,
    B = scenario_B,
    C = scenario_C,
    D = scenario_D
  ),
  uncertainty = TRUE,
  output = "summary",
  n_samples = 100,
  seed = 123
)

comp_pred <- comp_pred |>
  dplyr::mutate(
    comparison = paste(
      scenario1,
      scenario2,
      sep = " vs. "
    )
  )

output = "samples"

comp_pred_samples <- compare_predictions(
  fit = fit,
  profiles = list(
    A = scenario_A,
    B = scenario_B,
    C = scenario_C,
    D = scenario_D
  ),
  uncertainty = TRUE,
  output = "samples",
  n_samples = 100,
  seed = 123
)

comp_pred_samples <- comp_pred_samples |>
  dplyr::mutate(
    comparison = paste(
      scenario1,
      scenario2,
      sep = " vs. "
    )
  )

comparison_order <- comp_pred |>
  dplyr::arrange(diff) |>
  dplyr::pull(comparison)

comp_pred <- comp_pred |>
  dplyr::mutate(
    comparison = factor(
      comparison,
      levels = comparison_order
    )
  )

comp_pred_samples <- comp_pred_samples |>
  dplyr::mutate(
    comparison = factor(
      comparison,
      levels = comparison_order
    )
  )

sc_simu_plot <- ggplot(
  comp_pred,
  aes(
    x = comparison,
    y = diff
  )
) +
  geom_jitter(
    data = comp_pred_samples,
    aes(
      x = comparison,
      y = diff,
      color = sample
    ),
    alpha = 0.4,
    size = 2,
    width = 0.20,
    height = 0,
    inherit.aes = FALSE
  ) +
  scale_color_viridis_b() +
  geom_pointrange(
    aes(
      ymin = diff_lower,
      ymax = diff_upper
    ),
    color = "black",
    linewidth = 0.6
  ) +
  geom_hline(
    yintercept = 0,
    linetype = 2,
    color = "gray40",
    linewidth = 0.5
  ) +
  theme_bw() +
  labs(
    x = "Scenario comparison",
    y = "Difference in expected response"
  ) +
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    ),
    legend.position = "none",
    axis.text.x = element_text(
      angle = 45,
      hjust = 1
    )
  )

sc_simu_plot

scenario_low_rain <- list(
  tmean = exposure_tmean_A,
  rain = simulate_exposures(
    max_lag = 85,
    mode = "profile",
    seed = 501,
    background = list(
      dist = "ar1",
      mean = 2,
      sd = 1,
      phi = 0.50
    ),
    n = 1
  ),
  wetness = exposure_wetness_A
)

scenario_high_rain <- list(
  tmean = exposure_tmean_A,
  rain = simulate_exposures(
    max_lag = 85,
    mode = "profile",
    seed = 502,
    background = list(
      dist = "ar1",
      mean = 15,
      sd = 3,
      phi = 0.50
    ),
    n = 1
  ),
  wetness = exposure_wetness_A
)

rainfall_comparison <- compare_predictions(
  fit = fit,
  profiles = list(
    low_rain = scenario_low_rain,
    high_rain = scenario_high_rain
  ),
  uncertainty = TRUE,
  output = "summary",
  n_samples = 100,
  seed = 123
)

rainfall_comparison
#>   scenario1 scenario2     pred1   pred1_sd pred1_lower pred1_upper     pred2
#> 1  low_rain high_rain 0.8852602 0.04879873   0.7618158   0.9553559 0.9981504
#>       pred2_sd pred2_lower pred2_upper      diff    diff_sd diff_lower
#> 1 0.0009853875   0.9952408   0.9991543 0.1133656 0.04937617 0.04094392
#>   diff_upper  abs_diff abs_diff_sd abs_diff_lower abs_diff_upper
#> 1  0.2372621 0.1133656  0.04937617     0.04094392      0.2372621
#>   percentage_point_change percentage_point_change_sd
#> 1                11.33656                   4.937617
#>   percentage_point_change_lower percentage_point_change_upper
#> 1                      4.094392                      23.72621
#>   relative_defined_fraction relative_change_defined    ratio   ratio_sd
#> 1                         1                    TRUE 1.128128 0.06699667
#>   ratio_lower ratio_upper percent_change percent_change_sd percent_change_lower
#> 1    1.042867    1.311432       12.81284          6.699667             4.286725
#>   percent_change_upper
#> 1             31.14316

Interpreting simulation and prediction outputs

Simulation and prediction tools provide related but distinct outputs:

• simulate_exposures() creates one complete exposure profile;
• simulate_ranges() constructs combinations of candidate exposure values;
• simulate_scenarios() predicts outcomes over multiple scenarios;
• predict_outcomes() predicts the expected response for complete profiles;
• compare_exposures() compares the exposure profiles themselves;
• compare_predictions() compares their model-based expected outcomes.

Before interpreting a result, verify:

1. whether the exposure profiles are chronological;
2. whether every profile contains max_lag + 1 observations;
3. whether all exposures required by the fitted model were supplied;
4. whether scenario values are supported by the observed data;
5. whether multiple exposures changed simultaneously;
6. whether uncertainty was propagated;
7. whether the output contains summaries or individual samples;
8. whether the prediction is population-level or group-specific;
9. whether the interpretation is associational or supported by a causal design.

Summary

This section demonstrated how EpiExposure transforms a fitted DLNM into a framework for evaluating complete exposure profiles.

The workflow allows users to:

• generate chronological exposure profiles;
• add structured exposure patterns;
• predict expected disease outcomes;
• propagate parameter uncertainty;
• evaluate scenario grids;
• compare exposure profiles;
• compare predicted responses.

These tools extend inference from the fitted exposure-lag surface to explicit, model-based comparisons of complete environmental profiles.

Next steps

The next section translates predicted disease outcomes into practical outcomes such as yield and economic losses using the simulate_losses() function.

These tools connect exposure-lag-response relationships with agronomic and economic consequences.