Transmission Metrics

Calculate transmission metrics to estimate epidemic potential, velocity, and spread

Effective Reproduction

The Effective Reproduction Number (Re) is the average number of secondary premises infected by a source per day. Re is used to estimate epidemic potential.

Code

## exclude index cases (source == 0)
transmissions <- merge %>% 
  filter(source != 0)

## find the number of transmissions per source for each iteration and infect_day
iteration_Re <- transmissions %>%
  group_by(region, scenario_type, preclinical, iteration, infect_day, source) %>%
  summarize(num_transmissions = n(), .groups = "drop") 

iteration_mean_Re <- iteration_Re %>%
  group_by(region, scenario_type, preclinical, iteration, infect_day) %>%
  summarize(Re = mean(num_transmissions), .groups = "drop")

## summarize daily Re across iterations
daily_Re_summary <- iteration_mean_Re %>%
  group_by(region, scenario_type, preclinical, infect_day) %>%
  summarize(
    mean = mean(Re, na.rm = TRUE),
    q05 = quantile(Re, probs = 0.05, na.rm = TRUE),
    q25 = quantile(Re, probs = 0.25, na.rm = TRUE),
    q50 = quantile(Re, probs = 0.50, na.rm = TRUE),
    q75 = quantile(Re, probs = 0.75, na.rm = TRUE),
    q95 = quantile(Re, probs = 0.95, na.rm = TRUE),
    .groups   = "drop")

## Daily summary
daily_Re <- daily_Re_summary %>%
  group_by(region, scenario_type, preclinical) %>%
  complete(infect_day = 10:365, 
            fill = list(mean = 0, q05 = 0, 
                        q25 = 0, q50 = 0, q75 = 0, q95 = 0)) %>%
  ungroup()

Daily Summaries

Central U.S.
Eastern U.S.

View summary table of daily Re in the central region

Code

## Filter to central region
daily_Re_central <- daily_Re %>%
  filter(region == "central")

## Filter to scenario
daily_Re_central_select <- daily_Re_central %>%
  filter(scenario_type == "suboptimal") %>%
  filter(preclinical == "3")

## Select and order columns to display
daily_Re_central_select <- daily_Re_central_select[c("infect_day", "mean", "q05", "q25", "q50", "q75", "q95")]

head(daily_Re_central_select)

infect_day	mean	q05	q25	q50	q75	q95
10	0.000000	0	0	0	0	0.000
11	2.155172	1	1	2	3	4.150
12	1.771429	1	1	1	2	4.000
13	1.955095	1	1	1	2	5.000
14	1.945014	1	1	1	2	5.425
15	1.886877	1	1	1	2	5.000

View summary table of daily Re in the eastern region

Code

## Filter to eastern region
daily_Re_eastern <- daily_Re %>%
  filter(region == "eastern")

## Filter to scenario
daily_Re_eastern_select <- daily_Re_eastern %>%
  filter(scenario_type == "suboptimal") %>%
  filter(preclinical == "3")

## Select and order columns to display
daily_Re_eastern_select <- daily_Re_eastern_select[c("infect_day", "mean", "q05", "q25", "q50", "q75", "q95")]

head(daily_Re_eastern_select)

infect_day	mean	q05	q25	q50	q75	q95
10	0.000000	0	0	0	0.000	0.0
11	2.153846	1	1	2	3.000	4.0
12	1.971014	1	1	2	3.000	4.6
13	2.043243	1	1	1	3.000	5.0
14	1.878851	1	1	1	2.375	4.5
15	2.013787	1	1	1	3.000	5.0

Plot Effective Reproduction

The dashed line is at Re = 1

Central U.S.
Eastern U.S.

Code

## Prepare data
daily_Re <- daily_Re_central %>%
  mutate(
    preclinical = recode(preclinical,
                               "0" = "None",
                               "1" = "1 Day",
                               "2" = "2 Days",
                               "3" = "3 Days",
                               "6" = "6 Days")
  ) %>%
  mutate(
    scenario_type = recode(scenario_type,
                               "optimal" = "Optimal",
                               "suboptimal" = "Suboptimal",
                               "low-virulence" = "Low-Virulence")
  )  

daily_Re$preclinical <- ordered(factor(daily_Re$preclinical), 
                                      c("None", "1 Day", "2 Days", "3 Days", "6 Days"))

## Assign colors
myCols <- kelly()

## Plot data
ggplot(daily_Re, aes(infect_day, q50, group = preclinical, col = preclinical)) +
  geom_hline(yintercept = 1, col="black", linewidth=0.5, linetype = "dashed") +
  geom_smooth(method = "loess", span=0.05, se=FALSE) +
  scale_color_manual(values = c(
    "None" = myCols[3],
    "1 Day" = myCols[4],
    "2 Days" = myCols[6],
    "3 Days" = myCols[7],
    "6 Days" = myCols[10]
    )) +
  facet_wrap(~ scenario_type, ncol=1) +
  labs(
    x = "Simulation Day",
    y = "Median Effective Reproduction (Re)",
    col = "Preclinical Duration",
    title = " "
  ) +
  theme_minimal() +
  theme(
    plot.margin    = unit(c(0.25, 0.25, 0.25, 0.25), "cm"),
    legend.position = "bottom",
    strip.text     = element_text(size = 18, face = "bold", color = "gray40"),
    axis.title.x   = element_text(size = 22, face = "bold"),
    axis.title.y   = element_text(size = 20, face = "bold"),
    axis.text.x    = element_text(size = 18, face = "bold"),
    axis.text.y    = element_text(size = 18, face = "bold"),
    legend.key.width = unit(3, "line"),
    legend.key.height = unit(1, "line"),
    legend.text = element_text(size = 16, face = "bold"),
    legend.title = element_text(size = 18, face = "bold"),
    plot.title     = element_text(size = 22, face = "bold", hjust = 0.5)
  )

Code

## Prepare data
daily_Re <- daily_Re_eastern %>%
  mutate(
    preclinical = recode(preclinical,
                               "0" = "None",
                               "1" = "1 Day",
                               "2" = "2 Days",
                               "3" = "3 Days",
                               "6" = "6 Days")
  ) %>%
  mutate(
    scenario_type = recode(scenario_type,
                               "optimal" = "Optimal",
                               "suboptimal" = "Suboptimal",
                               "low-virulence" = "Low-Virulence")
  ) 

daily_Re$preclinical <- ordered(factor(daily_Re$preclinical), 
                                      c("None", "1 Day", "2 Days", "3 Days", "6 Days"))

## Assign colors
myCols <- kelly()

## Plot data
ggplot(daily_Re, aes(infect_day, q50, group = preclinical, col = preclinical)) +
  geom_hline(yintercept = 1, col="black", linewidth=0.5, linetype = "dashed") +
  geom_smooth(method = "loess", span=0.05, se=FALSE) +
  scale_color_manual(values = c(
    "None" = myCols[3],
    "1 Day" = myCols[4],
    "2 Days" = myCols[6],
    "3 Days" = myCols[7],
    "6 Days" = myCols[10]
    )) +
  facet_wrap(~ scenario_type, ncol=1) +
  labs(
    x = "Simulation Day",
    y = "Median Effective Reproduction (Re)",
    col = "Preclinical Duration",
    title = " "
  ) +
  theme_minimal() +
  theme(
    plot.margin    = unit(c(0.25, 0.25, 0.25, 0.25), "cm"),
    legend.position = "bottom",
    strip.text     = element_text(size = 18, face = "bold", color = "gray40"),
    axis.title.x   = element_text(size = 22, face = "bold"),
    axis.title.y   = element_text(size = 20, face = "bold"),
    axis.text.x    = element_text(size = 18, face = "bold"),
    axis.text.y    = element_text(size = 18, face = "bold"),
    legend.key.width = unit(3, "line"),
    legend.key.height = unit(1, "line"),
    legend.text = element_text(size = 16, face = "bold"),
    legend.title = element_text(size = 18, face = "bold"),
    plot.title     = element_text(size = 22, face = "bold", hjust = 0.5)
  )

Note

Example code demonstrates the workflow for the Eastern U.S. region. Central U.S. figures are included in the regional panels below. Select the tabs to view regional comparisons.

Epidemic Velocity

The epidemic velocity is the distance (km) of spread per day. The epidemic velocity is estimated by first calculating the daily spread distances.

Calculate Daily Distances

compile_daily_summary() calculates distances (km) between source farms and those infected for each iteration, then calculates the percentiles by day across all iterations to get the average statistics.

Code

daily_summary <- compile_daily_summary(region_paths, reference)

Compare Scenarios

Compare quantiles across all scenarios and prepare to calculate wave velocity.

Code

## Get quantiles
daily_all <- daily_summary$combined_summary

## Compare scenarios
all_summary_compare <- daily_all %>%
  mutate(
    scenario_type = str_to_title(word(scenario, 1, sep = "/")),
    preclinical = word(scenario, 2, sep = "/") %>% 
      str_replace_all("-days?-preclinical", "") %>% 
      str_trim() %>% 
      paste("Days Preclinical", .) %>% 
      str_to_title()
  )

Code

all_summary_compare <- all_summary_compare_central %>%
  mutate(
    preclinical = recode(preclinical,
                               "Days Preclinical 0" = "None",
                               "Days Preclinical 1" = "1 Day",
                               "Days Preclinical 2" = "2 Days",
                               "Days Preclinical 3" = "3 Days",
                               "Days Preclinical 6" = "6 Days")
  )

all_summary_compare$scenario_type <- ordered(factor(all_summary_compare$scenario_type), 
                                      c("Optimal", "Suboptimal", "Low-Virulence"))

all_summary_compare$preclinical <- ordered(factor(all_summary_compare$preclinical), 
                                      c("None", "1 Day", "2 Days", "3 Days", "6 Days"))

myCols <- kelly()

ggplot(all_summary_compare, aes(infect_day, log(q50), 
                               group = scenario, col = preclinical)) +
  geom_line(linewidth = 0.4) +
  ylim(0, 6) +
  scale_color_manual(values = c(
    "None" = myCols[3],
    "1 Day" = myCols[4],
    "2 Days" = myCols[6],
    "3 Days" = myCols[7],
    "6 Days" = myCols[10]
    )) +
  facet_wrap(~ scenario_type, ncol=1) +
  labs(
    x = "Simulation Day",
    y = "Distance (log km)",
    col = "Preclinical Infectious Duration",
    title = " "
  ) +
  theme_minimal() +
  theme(
    plot.margin    = unit(c(0.25, 0.25, 0.25, 0.25), "cm"),
    legend.position = "bottom",
    strip.text     = element_text(size = 18, face = "bold", color = "gray40"),
    axis.title.x   = element_text(size = 22, face = "bold"),
    axis.title.y   = element_text(size = 22, face = "bold"),
    axis.text.x    = element_text(size = 18, face = "bold"),
    axis.text.y    = element_text(size = 18, face = "bold"),
    legend.key.width = unit(3, "line"),
    legend.key.height = unit(1, "line"),
    legend.text = element_text(size = 16, face = "bold"),
    legend.title = element_text(size = 18, face = "bold"),
    plot.title     = element_text(size = 22, face = "bold", hjust = 0.5)
  )

Code

all_summary_compare <- all_summary_compare %>%
  mutate(
    preclinical = recode(preclinical,
                               "Days Preclinical 0" = "None",
                               "Days Preclinical 1" = "1 Day",
                               "Days Preclinical 2" = "2 Days",
                               "Days Preclinical 3" = "3 Days",
                               "Days Preclinical 6" = "6 Days")
  )

all_summary_compare$scenario_type <- ordered(factor(all_summary_compare$scenario_type), 
                                      c("Optimal", "Suboptimal", "Low-Virulence"))

all_summary_compare$preclinical <- ordered(factor(all_summary_compare$preclinical), 
                                      c("None", "1 Day", "2 Days", "3 Days", "6 Days"))

myCols <- kelly()

ggplot(all_summary_compare, aes(infect_day, log(q50), 
                               group = scenario, col = preclinical)) +
  geom_line(linewidth = 0.4) +
  ylim(0, 6) +
  scale_color_manual(values = c(
    "None" = myCols[3],
    "1 Day" = myCols[4],
    "2 Days" = myCols[6],
    "3 Days" = myCols[7],
    "6 Days" = myCols[10]
    )) +
  facet_wrap(~ scenario_type, ncol=1) +
  labs(
    x = "Simulation Day",
    y = "Distance (log km)",
    col = "Preclinical Infectious Duration",
    title = " "
  ) +
  theme_minimal() +
  theme(
    plot.margin    = unit(c(0.25, 0.25, 0.25, 0.25), "cm"),
    legend.position = "bottom",
    strip.text     = element_text(size = 18, face = "bold", color = "gray40"),
    axis.title.x   = element_text(size = 22, face = "bold"),
    axis.title.y   = element_text(size = 22, face = "bold"),
    axis.text.x    = element_text(size = 18, face = "bold"),
    axis.text.y    = element_text(size = 18, face = "bold"),
    legend.key.width = unit(3, "line"),
    legend.key.height = unit(1, "line"),
    legend.text = element_text(size = 16, face = "bold"),
    legend.title = element_text(size = 18, face = "bold"),
    plot.title     = element_text(size = 22, face = "bold", hjust = 0.5)
  )

Iteration Metrics

Compare scenarios by looking for patterns at the level of individual iterations. iteration_metrics() will use the other output in daily_summary to look at trends between optimal and suboptimal responses.

The function returns iteration specific metrics:

auc_log: The amount of area (total) under the plotted line (distance x day) above on a log scale, i.e., the Area Under Curve (AUC). This represents the total distance covered during the outbreak`s spread (on the log scale).
peak_spread: The maximum distance spread in any one day (log scale).
peak_day: The infect day that the maximum spread distance (peak_spread) occurred.

Code

## Calculate daily iteration metrics
iteration_metrics <- iteration_metrics(daily_summary$daily_distances) %>%
  select(-c(scenario))

## Capitalize scenario labels
iteration_metrics$scenario_type <- str_to_title(iteration_metrics$scenario_type)

## Order scenario levels
iteration_metrics$scenario_type <- ordered(factor(iteration_metrics$scenario_type), 
                                      c("Optimal","Suboptimal", "Low-virulence"))

View sample of iteration metrics

Code

## Filter data to scenario
iteration_metrics_select <- iteration_metrics %>%
  filter(scenario_type == "Suboptimal") %>%
  filter(preclinical == "3")

## Select and order columns to display
iteration_metrics_select <- iteration_metrics_select[c("iteration", "auc_log", "peak_spread", "peak_day")]

head(iteration_metrics_select)

iteration	auc_log	peak_spread	peak_day
1	57.56127	5.302539	18
2	231.03620	5.451650	25
3	901.49964	5.527616	89
4	162.08489	7.021728	36
5	244.14070	5.940245	12
6	197.69113	5.667478	56

Cumulative Spread

Drop the low-virulence scenarios to compare optimal and suboptimal detection scenarios. Plot the estimated Area Under the Curve (AUC). auc_log represents the total cumulative spread distance, for each scenario by iteration, on a log scale.

Code

## Drop the low-virulence scenario
iteration_metrics_no_low <- iteration_metrics %>%
  filter(scenario_type != "Low-Virulence")

## Plot AUC
ggplot(iteration_metrics_no_low, 
       aes(x = preclinical, y = auc_log, color = scenario_type)) +
  geom_point(shape=1, alpha = 0.3, size = 2) + # individual iterations
  geom_smooth(method = "lm", se = FALSE, linewidth=1.2) + # trend
  facet_grid(. ~ scenario_type) +
  scale_color_manual(values = c("Suboptimal" = "#74add1", "Optimal" = "orange2")) +
  labs(
    x = "Preclinical Infectious Duration (days)",
    y = "Cumulative Spread (log AUC)",
    title = " "
  ) +
  theme_minimal() +
  theme(
    plot.margin    = unit(c(0.25, 0.25, 0.25, 0.25), "cm"),
    legend.position = "none",
    strip.text     = element_text(size = 18, face = "bold", color = "gray40"),
    axis.title.x   = element_text(size = 22, face = "bold"),
    axis.title.y   = element_text(size = 22, face = "bold"),
    axis.text.x    = element_text(size = 18, face = "bold"),
    axis.text.y    = element_text(size = 18, face = "bold"),
    plot.title     = element_text(size = 22, face = "bold", hjust = 0.5)
  )

Significance Test

Fit a linear model with an interaction term to quantify the influence of preclinical infectious duration (preclinical) and detection scenario (scenario_type) on cumulative outbreak spread (auc_log).

Linear Model

Central U.S.
Eastern U.S.

Code

model_auc_central <- lm(auc_log ~ preclinical * scenario_type, data = iteration_metrics_no_low_central)

summary(model_auc_central)


Call:
lm(formula = auc_log ~ preclinical * scenario_type, data = iteration_metrics_no_low_central)

Residuals:
    Min      1Q  Median      3Q     Max 
-379.62  -92.16  -32.04   64.68  918.16 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                   79.564      5.467  14.554  < 2e-16 ***
preclinical                   58.007      2.753  21.070  < 2e-16 ***
scenario_type.L               55.472      7.731   7.175 8.85e-13 ***
preclinical:scenario_type.L   42.508      3.893  10.918  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 172.8 on 3391 degrees of freedom
  (301 observations deleted due to missingness)
Multiple R-squared:  0.294, Adjusted R-squared:  0.2934 
F-statistic: 470.8 on 3 and 3391 DF,  p-value: < 2.2e-16

Code

model_auc <- lm(auc_log ~ preclinical * scenario_type, data = iteration_metrics_no_low)

summary(model_auc)


Call:
lm(formula = auc_log ~ preclinical * scenario_type, data = iteration_metrics_no_low)

Residuals:
    Min      1Q  Median      3Q     Max 
-316.97  -79.18  -28.81   47.52  855.07 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                   78.708      4.674  16.840  < 2e-16 ***
preclinical                   45.739      2.373  19.276  < 2e-16 ***
scenario_type.L               48.654      6.610   7.361 2.28e-13 ***
preclinical:scenario_type.L   32.891      3.356   9.802  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 150.3 on 3429 degrees of freedom
  (280 observations deleted due to missingness)
Multiple R-squared:  0.2637,    Adjusted R-squared:  0.263 
F-statistic: 409.3 on 3 and 3429 DF,  p-value: < 2.2e-16

ANOVA

Central U.S.
Eastern U.S.

Code

anova(model_auc_central)

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
preclinical	1	12687728	12687727.73	424.8123	0
scenario_type	1	25933670	25933669.92	868.3148	0
preclinical:scenario_type	1	3560019	3560019.23	119.1971	0
Residuals	3391	101277875	29866.67	NA	NA

Code

anova(model_auc)

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
preclinical	1	8014986	8014986.44	354.99578	0
scenario_type	1	17537524	17537523.84	776.76326	0
preclinical:scenario_type	1	2169123	2169123.39	96.07373	0
Residuals	3429	77418915	22577.69	NA	NA

Peak Spread

As with cumulative spread above, comparisons can be made using the peak spread rates.

Code

ggplot(iteration_metrics_no_low, aes(x = preclinical, y = peak_day, color = scenario_type)) +
  geom_point(shape=1, alpha = 0.3, size = 2) + # individual iterations
  geom_smooth(method = "lm", se = FALSE, linewidth=1.2) +
  facet_grid(. ~ scenario_type) +
  scale_color_manual(values = c("Suboptimal" = "#74add1", "Optimal" = "orange2")) +
  labs(
    x = "Preclinical Infectious Duration (days)",
    y = "Day of Peak Transmission",
    title = " "
  ) +
  theme_minimal() +
  theme(
    plot.margin    = unit(c(0.25, 0.25, 0.25, 0.25), "cm"),
    legend.position = "none",
    strip.text     = element_text(size = 18, face = "bold", color = "gray40"),
    axis.title.x   = element_text(size = 22, face = "bold"),
    axis.title.y   = element_text(size = 22, face = "bold"),
    axis.text.x    = element_text(size = 18, face = "bold"),
    axis.text.y    = element_text(size = 18, face = "bold"),
    plot.title     = element_text(size = 22, face = "bold", hjust = 0.5)
  )

Significance Test

Fit a linear model with an interaction term to quantify the influence of preclinical and scenario_type on the maximum spread rate and it’s (peak_day).

Linear Model

Code

model_peak <- lm(peak_day ~ preclinical * scenario_type, data = iteration_metrics_no_low)

summary(model_peak)


Call:
lm(formula = peak_day ~ preclinical * scenario_type, data = iteration_metrics_no_low)

Residuals:
    Min      1Q  Median      3Q     Max 
-43.909 -13.940  -6.919   5.852 287.081 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  26.3329     0.8619  30.553  < 2e-16 ***
preclinical                   4.7061     0.4506  10.444  < 2e-16 ***
scenario_type.L               6.5160     1.2189   5.346 9.54e-08 ***
preclinical:scenario_type.L   4.6434     0.6373   7.286 3.87e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 30.38 on 3709 degrees of freedom
Multiple R-squared:  0.1263,    Adjusted R-squared:  0.1256 
F-statistic: 178.7 on 3 and 3709 DF,  p-value: < 2.2e-16

ANOVA

Code

anova(model_peak)

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
preclinical	1	97130.83	97130.8281	105.25678	0
scenario_type	1	348554.42	348554.4173	377.71442	0
preclinical:scenario_type	1	48993.79	48993.7929	53.09261	0
Residuals	3709	3422660.79	922.7988	NA	NA