submission_2.Rmd

title: "ASF Challenge - submission 2"
author: "Cirad-Inrae team"
date: "`r format(Sys.Date(), '%e %B, %Y')`"
# knit: (function(inputFile, encoding, knit_root_dir) {
#     rmarkdown::render(
#       inputFile,
#       encoding = encoding,
#       knit_root_dir = knit_root_dir,
#       output_dir = "../report",
#       output_format = "pdf_document"
#     )
#   })
output:  bookdown::pdf_document2
classoption: "a4paper"
bibliography: asf_challenge.bib

source("src/packages.R")
source("src/functions.R")
message("Assuming interactive session")

knitr::opts_chunk$set(
  echo = FALSE,
  cache = FALSE
)

theme_set(theme_bw())

## Explicitly name all targets we want access to
loadd(
  admin,
  animal_movements,
  # data_structure,
  # distMatPigs,
  # farm_exposition_D80,
  farm_pred_inf_ship_CSIF_D50,
  farm_pred_inf_ship_CSIF_D80,
  fence,
  fm_farm_inf_D80,  # model prob farm infection
  # herd_cullings_CSIF_D80,
  # herd_cullings_CSPZ_D80,
  # herd_cullings_CSWB_D80,
  # herd_cullings_CSTR_D80,
  # hexAll,
  huntZone,
  # huntedWBisland,
  huntedWBzoom,
  IWB_scenarios,
  # landcover,
  # lcMap,
  # lc_lite,
  # lc_lite_r,
  mcmcFileModel10,
  # movementAnimation,
  # movementAnimationHack,
  # movementFileName,
  # negative_hunt_samples,
  # obs_days_D80,         ## Number of observation days for the exposure process
  # outbreaks,
  # outbreaks_at_D50,
  outbreaks_at_D80,
  pig_sites,
  pig_sites_pred_CSIF_D50,
  pig_sites_pred_CSIF_D80,
  pig_sites_risks_D80,
  pigSitesExposure,     # Based on IWB simulations up to day 140
  tab_farms_pred_CSIF_D50,
  tab_farms_pred_CSIF_D80,
  tab_farms_pred_CSPZ_D80,
  tab_farms_pred_CSWB_D80,
  tab_farms_pred_CSTR_D80,
  scenariosWBR2,
  map_expectedI_scenario1_140,
  map_expectedI_scenario2_140,
  map_expectedI_scenario3_140,
  map_expectedI_scenario4_140,
  tab_probAS,
  # tmAdminPigTypesAndOutbreaks,
  # tmProportionAgricultureMap,
  # tmProportionForestCoverMap,
  # tblMoves,
  tmAdmin,
  # tmRowAdmin,
  wb_case_density_D80,
  # wildBoarObsAll,
)

Introduction

In this second report we analyse data made available on day 80 -- three months since the ASF epidemic started on Merry Island. To control spread among wild boar (WB), a fence was erected on day 60 and it was expected that it would isolate the infected WB population. In parallel, an increasing hunting pressure was implemented in a buffer zone surrounding and including the fence area, to decrease wild boar densities and slow the spread of the virus. In addition, the inclusion of the army into the control efforts has increased the passive search detection rate. Many new ASF cases among WB have been detected, both through the testing of shot WB and through passive detection. During this time several new measures have been implemented, including:

1. An increased hunting pressure within a buffer zone of 15 km around the fenced area.
2. Surveillance and protection zones have been created around infected domestic pig farms and the shipment of animals out of these zones has been banned.
3. All domestic pigs on infected farms now have to be culled after the detection of a case.
4. Since the disease continues to spread, authorities are considering to apply new culling policies.
5. There is currently a question as to whether the size of the active search buffer around carcasses should be increased from 1 km to 2 km.

In the present report, we update our previous predictions and quantify the effectiveness of the following proposed strategies on the spread of the disease:

1. Presence of fences, with or without an increased hunting pressure in the fenced area and newly created 15km buffer zone.
2. Increasing the active search radius around WB carcasses from 1km to 2km.
3. Culling all pig herds in protection zones.
4. Culling pig herds located at 3 km of WB carcasses.
5. Culling all pig herds that have traded pigs with infected farms less than 3 weeks before detection.

Assesment of previous predictions

The predictions of our previous report are contrasted with subsequent observations in figure @ref(fig:farms-pred-vs-realised).

outbreaks_predD50_vs_realised <-
  full_join(
    outbreaks_at_D80 %>%
      filter(HOST == "pig herd") %>%
      st_drop_geometry() %>%
      mutate(Infection = ifelse(DATE.SUSP <= 50, "Initial", "New"))
    ,
    tab_farms_pred_CSIF_D50 %>%
      select(-size) %>%
      mutate(predicted = p_pred > 0.05) %>%
      filter(predicted),
    by = c(ID = "id")
  ) %>%
  replace_na(
    list(
      Infection = "Not infected",
      predicted = FALSE
    )
  ) %>%
  mutate(Infection = fct_inorder(Infection)) %>%
  left_join(
    pig_sites,
    by = c(ID = "population_id")
  ) %>%
  st_sf()

cap <- "Predicted risks made at day 50 (using updated methods) and observed infections at pig sites up to day 80. Points outlined in red were predicted with probability lower than 5%. The square (indoor) farm was not predicted at risk with the previous methods."

tm_shape(
  admin,
  bbox = st_bbox(
    outbreaks_predD50_vs_realised %>%
      st_buffer(dist = 5e4)
  )
) +
  tm_layout(bg.color = "lightblue") +
  # tm_grid(col = "grey", labels.size = 1) +
  tm_scale_bar(
    breaks = c(0, 50, 100, 150),
    text.size = 0.9,
    position = c("right", "bottom")
  ) +
  tm_polygons() +
  ## Pig sites
  tm_shape(
    outbreaks_predD50_vs_realised %>% filter(predicted)
  ) +
  tm_bubbles(size = "p_pred", scale = 1, col = "Infection") +
  tm_shape(
    outbreaks_predD50_vs_realised %>% filter(!predicted)
  ) +
  tm_bubbles(
    border.col = "darkred", col = "Infection", size = .1, shape = "practice",
    border.lwd = 2, legend.col.show = FALSE
  ) +
  ## Animal shipments at risk
  # tm_shape(pig_sites_pred_CSIF_D50) +
  #   tm_dots(col = "FarmType", size=0.3, palette = "viridis") +
  tm_shape(
    animal_movements %>%
      filter(date >= 50 - detection_delay, date <= 50) %>%
      filter(
        source %in%
          filter(pig_sites_pred_CSIF_D50, susceptible_wb, !detected)$population_id,
        dest %in% farm_pred_inf_ship_CSIF_D50$dest
      )
  ) +
  tm_lines() +
  ## Fences
  tm_shape(c(st_geometry(fence), st_geometry(huntZone))) +
  tm_borders()

unpredicted_infections <-
  outbreaks_predD50_vs_realised %>%
  filter(!predicted, Infection == "New", is_outdoor == 0)

unpredicted_infections %>%
  st_drop_geometry() %>%
  select(ID, DET, starts_with("DATE"), FarmType, practice, size) %>%
  kable(
    format = "latex",
    booktabs = TRUE,
    caption = "List of newly infected sites not identified as at risk."
  )

susp_contact_case <-
  filter(pig_sites_pred_CSIF_D50, susceptible_wb)[
  st_nearest_feature(
    unpredicted_infections,
    pig_sites_pred_CSIF_D50 %>% filter(susceptible_wb)
  ),
]

dist_to_susp_contact_case <-
  st_distance(
    unpredicted_infections,
    susp_contact_case
  ) %>% as.numeric

pig_sites_pred_CSIF_D80 %>%
  inner_join(
    susp_contact_case %>%
      st_drop_geometry() %>%
      select(population_id),
    by = "population_id"
  ) %>%
  st_drop_geometry() %>%
  select(Id = population_id, FarmType, practice, size, "day symptoms" = date_susp) %>%
  kable(
    format = "latex",
    booktabs = TRUE,
    caption = "Suspected contact case."
  )

There has been r nrow(unpredicted_infections) newly infected sites that we failed to identify as being at risk (Table @ref(tab:unpredicted-sites)). We had not identified this site because it was an in-door site, and thus we had excluded it in our previous analysis. However, it is located in close proximity (~r round(dist_to_susp_contact_case, -1) m) to another infected, albeit very small, outdoor finishing farm (Table @ref(tab:susp-contact-case)). This suggests a fomite transmission pathway, related to proximity to infected farms. Thus, in this report, we extend our risk analysis to include this fomites transmission pathway (section @ref(sectionPigFarms)).

In our first report, we had predicted that ASF would spread through areas of high wild boar density as a travelling wave (figure @ref(fig:predictedIWB)). Moreover, we had indicated that this travelling wave would likely:

1. be somewhat slow and broken in the west due to low wild boar densities;
2. spread rapidly to the east and north-east, although probably not beyond the fence;
3. pass the proposed fence in the south prior to day 60, from where further spread among large wild boar populations was likely.

We had concluded that it was this third possibility that pressented the greatest threat to pig production chains on Merry Island.

cap <- "Our prediction from report 1 concerning the density of infectious wild boar (DIWB) for days 60 (left) and 80 (right).
We had predicted that the wave front would remain within the fenced area in the north and east, but that it would pass the proposed fence in the south before day 60. "
knitr::include_graphics(here("figures/report1/infectiousBest_fence100_tFence60_t60-80_half.jpg"))

Our predictions proved to be acurate, as indicated by the distribution of ASF-positve cases among all hunted wild boar reported up to day 80 (figure @ref(fig:observedIWB)).

cap <- "Locations of all hunted and tested wild boar (points).
The fence (built day 60) is marked purple.
The red rectangle is the area of increased hunting pressure.
The shading of the administrative areas (polygons) corresponds to the 2019 hunting bag."
print(huntedWBzoom)

Given that the latest data support the model's predictive capactiy, we chose to maintain it's basic structure. However, we have made a few modifications to improve: it's fit to the data; the convergence of our MCMC algorithm; and to account for the increased hunting pressure around the fenced zone.

In this report we present our work in two sections: the first is an analysis of risk among pig farms; the second an analysis of disease spread among wild boar. We finish the report with a brief section outlining our main conclusions.

\clearpage

Pig farms {#sectionPigFarms}

Methods

The methods used for farm risk assessments were described in the previously submitted report with some modifications and additions described here.

In addition to the infection pathways considered in the last submission, we included fomites as a possible infection pathway, for instance due to veterniarians or visitors that may introduce the virus, via their equipment or vehicles, after visiting an infected farm.

We consider the risk higher when infected farms, or farms with high probability of being infected, are in close proximity to the focal site. We thus model this risk as a exponentially decreasing function of the distance between farms.

Since infected farms in the vicinity may not yet have been detected, we consider all neighbouring farms weighted by their exposition to infectious wild boar.

Another change from the previous submission is a correction of the estimation of the probability of infection via animal shipment. Previously we considered all the recorded shipments as potential transmission events. Now, we consider only shipments that have taken place in the last 15 days. Earlier shipments can't be transmission events, otherwise any subsequent infection at a destination farm would have been detected.

We now also account for the culling strategy of infected herds, which was disregarded in our first report since there were very few. Now we consider that culled herds can only be infected again after repopulation, which is 50 days after infection. This lowers the risk of infection due to wild boar exposure for some farms, which also impacts the risk of transmission via fomites.

Finally, we address the three proposed culling strategies and assess their efficacy based on their predicted impact. Specifically, we compute the expected number of infected pigs in the period from D80 to D110 by summing the size of the farms weighted by their estimated probability of infection through any of the considered pathways.

sites_area <-  pig_sites_pred_CSIF_D80 %>% filter(expo_inf_wb > 0)  # not only outdoor!! and also near other infected farms.

# sites_area <- pig_sites %>%
#       add_column(
#         expo = farm_exposition_D80
#       ) %>%
#       ## Include variable "detected"
#       left_join(
#         outbreaks_at_D80 %>%
#           st_drop_geometry() %>%
#           filter(HOST == "pig herd") %>%
#           select(population_id = ID) %>%
#           add_column(detected = TRUE),
#         by = "population_id"
#       ) %>%
#       replace_na(list(detected = FALSE)) %>%
#       ## Filter exposed sites and transform detected status into 0/1
#       ## We are considering only outdoor farms. Perhaps the type B, BF, F
#       ## also plays a role? Difficult to assess at this point.
#       filter(expo > 0)

distance_to_nearest_infected_site <-
  st_distance(
    pig_sites_risks_D80 %>%
      filter(!detected, susceptible_wb),
    pig_sites_risks_D80 %>%
      filter(detected == 1, susceptible_wb),
  ) %>%
  apply(1, min)

## Assume a exponentially decreasing exposition kernel
## with mean of 700 m (the only suspected case) and multiply
## by the probability of wild boar infection, taking 1
## for already infected farms.

## Calibrate the kernel so that the risk halves at 700m


pig_sites_pred_CSIF_D80 %>%
  mutate(
    p_fomites = prob_fomites(p_inf, unclass(st_distance(.)))
  ) %>%
  filter(p_fomites > 0.005) %>%
  ggplot(aes(p_inf, p_fomites)) +
  geom_point() +
  geom_point(
    data = function(x) filter(x, detected),
    colour = "orange"
  ) +
  geom_text_repel(
    data = function(x) filter(x, detected),
    aes(label = population_id),
    colour = "orange"
  ) +
  # scale_x_log10() +
  # scale_y_log10() +
  coord_fixed()

idx <-
  pig_sites_pred_CSIF_D80 %>%
  mutate(
    p_fomites = prob_fomites(p_inf, unclass(st_distance(.)))
  ) %>%
  slice_max(p_fomites, n = 5)

\clearpage

Results

Farm circulation under current strategy

cap <- "Region at risk for outdoor farm contamination due to proximity to infected wild boar.
Background colour represents a kernel density smoothing of wild boar case density.
Smaller red dots are the locations of WB cases. Bigger white
points are outdoor farms. Those already infected and detected by day 80
are marked with a black dot. kde indicates kernel density estimate,
and is used to quantify the exposition level of outdoor farms."

infected_sites_region <-
  outbreaks_at_D80 %>%
  filter(HOST == "pig herd") %>%
  select(ID) %>%
  st_set_agr("identity") %>%
  st_intersection(st_union(wb_case_density_D80))

ggplot() +
  ## WB density surface
  geom_sf(
    data = wb_case_density_D80,
    aes(colour = kde, fill = kde),
    lwd = 0
  ) +
  ## WB cases
  geom_sf(
    data = outbreaks_at_D80 %>% filter(HOST == "wild boar"),
    col = "red"
  ) +
  ## Farm sites (outdoor) predicted at risk
  geom_sf(
    data = pig_sites_pred_CSIF_D80 %>% filter(susceptible_wb),
    col = "gray90",
    size = 3
  ) +
  # ## Farm sites predicted at risk via fomites
  # geom_sf(
  #   data = pig_sites_pred_CSIF_D80 %>% filter(population_id %in% idx),
  #   col = "green",
  #   size = 3
  # ) +
  ## Infected farm sites
  geom_sf(
    data = infected_sites_region,
  ) +
  geom_text_repel(
    # Awesome trick to use ggrepel with geom_sf.
    # https://github.com/slowkow/ggrepel/issues/111#issuecomment-416853013
    data = infected_sites_region,
    aes(label = ID, geometry = geometry),
    stat = "sf_coordinates"
    # nudge_y = 1e4
  ) +
  ## Scales and aesthetics
  scale_fill_viridis_c() +
  scale_colour_viridis_c() +
  ggsn::scalebar(
    x.min = min(st_coordinates(wb_case_density_D80)[, "X"]),
    x.max = max(st_coordinates(wb_case_density_D80)[, "X"]),
    y.min = min(st_coordinates(wb_case_density_D80)[, "Y"]),
    y.max = max(st_coordinates(wb_case_density_D80)[, "Y"]),
    dist = 10,
    dist_unit = "km",
    transform = FALSE,
    # anchor = c(x = 860e3, y = 2e5),
    # nudge_y = -1e3,
    # height = .04,
    # border.size = .5,
    st.size = 3
  ) +
  labs(x = NULL, y = NULL)

cap <- "Probability of farm infection and detection (curve), estimated from pressence/absence data (points),
for outdoor farms in the region at risk, as a function of the exposition level."
pig_sites_risks_D80 %>%
  filter(susceptible_wb) %>%
  mutate(
    detected = as.numeric(detected)
  ) %>%
  left_join(
    outbreaks_at_D80 %>%
      filter(HOST == "pig herd") %>%
      st_drop_geometry() %>%
      transmute(
        population_id = ID,
        Infection = ifelse(DATE.SUSP <= 50, "Initial", "New")
      ),
    by = "population_id"
  ) %>%
  replace_na(list(Infection = "Not infected")) %>%
  ggplot(aes(expo_inf_wb, detected)) +
  geom_point(aes(group = population_id, colour = Infection)) +
  geom_line(
    data = data.frame(
      expo_inf_wb = seq(
        min(pig_sites_pred_CSIF_D80$expo_inf_wb),
        max(pig_sites_pred_CSIF_D80$expo_inf_wb),
        length = 101
      )
    ) %>%
      mutate(
        pred = predict(fm_farm_inf_D80, type = "response", newdata = .)
      ),
    aes(expo_inf_wb, pred)
  ) +
  labs(x = "Exposition level", y = "Detection probability") +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.x = element_blank()
  )

cap <- "Probability of infection in a number of days, given the daily infection probability."
expand.grid(
  p_d = seq(0.001, .050, length = 50),
  n_d = seq(45, 75)
) %>%
  mutate(
    p_p = p_period(p_d, n_d)
  ) %>%
  ggplot(aes(p_d, n_d)) +
  geom_tile(aes(fill = p_p)) +
  geom_contour(aes(z = p_p), breaks = .5, col = "white") +
  scale_fill_viridis_c() +
  labs(x = "Prob. infection in a day", y = "Number of days")

cap <- "Probability of detection due to exposition to infected wild boar during the prediction period (D80 - D110), for farms in the area at risk."

plot_area <- st_bbox(
  c(
    pig_sites_pred_CSIF_D80 %>% filter(susceptible_wb) %>% st_geometry(),
    pig_sites_pred_CSIF_D80 %>% filter(susceptible_wb) %>% st_geometry() %>% `+`(c(5e4, 0))
  ) %>%
    st_union() %>%
    st_buffer(dist = 1e4)
)

tm_shape(admin, bbox = plot_area) +
  tm_layout(bg.color = "lightblue") +
  # tm_grid(col = "grey", labels.size = 1) +
  tm_scale_bar(
    breaks = c(0, 10, 30, 60),
    text.size = 0.9,
    position = c("left", "top")) +
  tm_polygons() +
  tm_shape(pig_sites_pred_CSIF_D80 %>% filter(susceptible_wb)) +
  tm_bubbles(
    title.size = "Prob. of detection",
    size = "p_pred",
    col = "detected",
    scale = 2
  ) +
  tm_shape(fence) +
  tm_borders() +
  tm_shape(huntZone) +
  tm_borders() +
  tm_legend(
    bg.color = "white",
    bg.alpha = .7,
    frame = TRUE
  )

cap <- "Cartographic representation of observed animal shipments (at any time) with risk of ASFV transmission (from farm sites at risk)."
tmAdmin +
  tm_shape(pig_sites_pred_CSIF_D80 %>% filter(susceptible_wb)) +
  tm_bubbles(
    title.col = "Detected",
    col = "detected",
    size = 0.2,
    palette = "viridis"
  ) +
  tm_shape(
    pig_sites %>%
      right_join(
        farm_pred_inf_ship_CSIF_D80,
        by = c(population_id = "dest")
      )
  ) +
  tm_bubbles(
    title.size = "Exp.# of infectious shipments.",
    size = "e_inf_shipments"
  ) +
  tm_shape(
    animal_movements %>%
      filter(
        source %in%
          filter(pig_sites_pred_CSIF_D80, susceptible_wb)$population_id,
        dest %in% farm_pred_inf_ship_CSIF_D80$dest)
  ) +
  tm_lines() +
  tm_shape(fence) +
  tm_borders() +
  tm_legend(
    bg.color = "white",
    bg.alpha = .7,
    frame = TRUE
  )

cap <- "Probability of detection due to exposition to fomites during the prediction period (D80 - D110). The rectangles represent the fenced (inner) and the buffer zones (outer)."

sites_fomites <-
  pig_sites_pred_CSIF_D80 %>%
  mutate(
    p_fomites = prob_fomites(p_inf, unclass(st_distance(.)))
  ) %>%
  filter(p_fomites > 0.01 | detected)

plot_area <- st_bbox(
  c(
    sites_fomites %>% st_geometry(),
    sites_fomites %>% st_geometry() %>% `+`(c(7e4, 0))
  ) %>%
    st_union() %>%
    st_buffer(dist = 1e4)
)

tm_shape(admin, bbox = plot_area) +
  tm_layout(bg.color = "lightblue") +
  # tm_grid(col = "grey", labels.size = 1) +
  tm_scale_bar(
    breaks = c(0, 10, 30, 60),
    text.size = 0.9,
    position = c("left", "top")) +
  tm_polygons() +
  tm_shape(sites_fomites) +
  tm_bubbles(
    title.size = "Prob. of detection",
    size = "p_fomites",
    col = "detected",
    scale = 2
  ) +
  tm_shape(fence) +
  tm_borders() +
  tm_shape(huntZone) +
  tm_borders() +
  tm_legend(
    bg.color = "white",
    bg.alpha = .7,
    frame = TRUE
  ) +
  tm_layout(
    legend.position = c("right", "top")
  )

kable(
  tab_farms_pred_CSIF_D80 %>%
    filter(p_pred > 0.05, !detected) %>%
    select(-detected),
  caption = "Probability of infection in undetected farm sites with non-negligible probability (P > 0.05) in the next period.",
  format = "latex",
  digits = 3,
  booktabs = TRUE
)

cap <- "Probability of infection and detection by any possible route during the prediction period (D80 - D110). The rectangles represent the fenced (inner) and the buffer zones (outer)."

sites_all <-
  pig_sites_pred_CSIF_D80 %>%
  select(-detected, -size, -p_pred) %>%
  inner_join(
    bind_rows(
      tab_farms_pred_CSIF_D80 %>%
        filter(p_pred > 0.01, !detected),
      tab_farms_pred_CSIF_D80 %>%  # Missing the detection in the north
        filter(detected)
    ),
  by = c(population_id = "id")
)

plot_area <- st_bbox(
  c(
    sites_all %>% st_geometry(),
    sites_all %>% st_geometry() %>% `+`(c(7e4, 0))
  ) %>%
    st_union() %>%
    st_buffer(dist = 1e4)
)

tm_shape(admin, bbox = plot_area) +
  tm_layout(bg.color = "lightblue") +
  # tm_grid(col = "grey", labels.size = 1) +
  tm_scale_bar(
    breaks = c(0, 10, 30, 60),
    text.size = 0.9,
    position = c("right", "bottom")) +
  tm_polygons() +
  tm_shape(sites_all) +
  tm_bubbles(
    title.size = "Prob. of detection",
    size = "p_pred",
    col = "detected",
    scale = 2
  ) +
  tm_shape(fence) +
  tm_borders() +
  tm_shape(huntZone) +
  tm_borders() +
  tm_legend(
    bg.color = "white",
    bg.alpha = .7,
    frame = TRUE
  ) +
  tm_layout(
    legend.position = c("right", "top")
  )

Impact of additional culling measures

Table @ref(tab:exp-pigs-cs) shows our preliminary calculations of the impact we expect from each of the proposed four control strategies. The current strategy is culling all infected farms (IF). The alternative strategies\footnote{PZ: culling all pig herds in protection zones; WB: culling all pig herds located at < 3 km from positive wild boar carcasses; TR: culling all pig herds that have traded pigs with an infected farm less than 3 weeks before detection} are in addition to IF and are therefore at least as effective.

(exp_pigs_cs <-
  bind_rows(
    IF = tab_farms_pred_CSIF_D80,
    PZ = tab_farms_pred_CSPZ_D80,
    WB = tab_farms_pred_CSWB_D80,
    TR = tab_farms_pred_CSTR_D80,
    .id = "CS"
  ) %>%
  filter(!detected) %>%
  group_by(CS) %>%
  summarise(
    `Exp pigs` = sum(p_pred * size),
    .groups = "drop"
  )) %>%
  kable(
    caption = "Expected number of infected pigs between D80 and D110 under alternative control strategies. IF: Infected farms; PZ: Protection Zone; TR: Trade; WB: Wild boar.",
    format = "latex",
    digits = 0,
    booktabs = TRUE
  )

\clearpage

Wild boar

Methods

cap <- "Compartmental model for ASF spread in wild boar described in report 1."
knitr::include_graphics(here("figures/modelWBflowDiagram.jpg"))

The spread of ASF among wild boar was analysed with a modified version of the stochastic compartmental model first presented in report 1 (figure @ref(fig:wbModel)). The table of parameters and associated priors of the modified model are as follows.

Notation	Parameter	Prior
t_\text{Intro}	ASF introduction time	-(Beta(7,7) x 100)
x_\text{Intro}	ASF introduction (easting)	Unif(extent of island)
y_\text{Intro}	ASF introduction (northing)	Unif(extent of island)
\beta	Transmission rate	Exp(\lambda=0.1)
p_\text{Home}	Connectivity (prob. WB in home pixel)	Beta(10,10)
p_\text{AttractI}	Attractivity of living IWB relative to carcass	Beta(2,2)
\tau_\text{P}	Passive search detection rate	LogNorm(lmean=0, lsd=3)
\tau_\text{A}	Baseline active search detection rate	Exp(\lambda=10^{-7})
\tau_\text{Phz}	Augmentation of \tau_\text{P} in hunting zone	Exp(\lambda=10^{-7})
\tau_\text{Ahz}	Augmentation of \tau_\text{A} in hunting zone	Exp(\lambda=10^{-7})
p_\text{HuntY}	Prob. hunted in 1 year	Beta(330,330)
\omega_\text{Fence}	Efficacy of fence	Beta(2,2)

Note, semi-informative priors were choosen for \beta and p_\text{Home} to improve MCMC performance for these two very highly correlated parameters.

Improved active search

The original model neglected the localised nature of active searches, and had assumed that the active search detection rate was uniform in space and time. In our modified model an active search affects only:

1. the pixel within which an infected wild boar was detected -- we call this the central pixel;
2. the pixels immediately adjacent to the central pixel -- i.e. neighbouring pixels.

Recall that our model uses a hexagonal grid with a 5km distance between pixel centroids. Thus, we utilise approximations at this resolution to investigate the difference between 1km and 2km active search radii.

Note that \tau_A is the initial value given to the active search detection rate within a central pixel when using a 1km search radius. We further assume:

1. a 7 day delay between a case confirmation date and active search initialisation -- this delay is consistently observed in the data;
2. once initialised, active search detection rates decrease each day by a multiplicative factor of (1-\frac{1}{7});
3. for a 1km search radius, the daily detection probability throughout a central pixel is initialised at p_{C1} = 1 - \exp(-\tau_A);
4. for a 1km search radius, the daily detection probability throughout a neighbouring pixel is initialised as p_{N1} = \frac{\omega_{N1}}{\omega_{C1}} p_{C1};
5. for a 2km search radius, the daily detection probability throughout central pixels and neighbouring pixels are p_{C2} = \frac{\omega_{C2}}{\omega_{C1}} p_{C1} and p_{N2} = \frac{\omega_{N2}}{\omega_{C1}} p_{C1} respectively;
6. the weights \omega_{Cx} and \omega_{Nx} indicate the expected proportions of a circle of radius x to fall within the central pixel, or a given neighbouring pixel, respectively;
7. when cases are found in multiple neighbours, detection rates accumulate in a pixel up to a maximum value of \tau_{Ax}=-\log(1-p_{Cx}).

Note, the above procedure is identical inside the zone of increased hunting pressure, except that the initial detection rate is augmented to \tau_A + \tau_{Ahz} for a central pixel under the 1km search radius scenario.

MCMC and simulation

The analysis was peformed in two steps.

1. First, we perfrom a parameter estimation step, fitting the model to the newly available data via Markov chain Monte Carlo (MCMC).
2. Secondly, we perfrom a prediction step to explore the relative efficacy of eight alternative control scenarios.

In the parameter estimation step, we employed the same Monte-Carlo likelihoods detailed in report 1. i.e. for each likelihood calculation, 500 simulations were used to estimate the expected number of cases for each of four observation types (hunt -ve, hunt +ve, active search and passive search) in each pixel, aggregated into ten-day bins. A series of 'burn-in' runs were used to tune the algorithm. Thereafter, inference was based on the output of 5 different MCMCs, each providing 4000 samples.

In the prediction step, we considered eight control scenarios (Table @ref(tab:iwbControlScenarios)) that cover various combinations of using / not using: the fence; an increase in hunting and passive search rate; a 1km or 2km active search radius. For each control scenario 10000 simulations were performed -- each based on a parameter set randomly selected from the 20000 lines of MCMC output. Each simulation ran from the estimated disease introduction date to day 140, thus providing a 40 day projection beyond day 80.

kable(scenariosWBR2,
    caption = "Table of the eight tested scenarios for controlling infectious wild boar",
    format = "latex",
    digits = 3,
    booktabs = TRUE
  )

Results

The expected density of infectious wild boar (IWB) under control scenarios 1-4 is pressented in figures @ref(fig:predicted-IWB-day140-scenario1), @ref(fig:predicted-IWB-day140-scenario2), @ref(fig:predicted-IWB-day140-scenario3) and @ref(fig:predicted-IWB-day140-scenario4). Note: the lightest yellow colour indicates areas where the expected IWB density is less than 1, i.e. areas where only a minority of the stochastic simulations indicated the pressence of IWB.

cap <- "Expected densities of infectious wild boar on day 140 under control scenario 1 (current control measures): fence (purple), increaed hunting zone (brown) and 1km active search (not shown)."
print(map_expectedI_scenario1_140)

Figure @ref(fig:predicted-IWB-day140-scenario1) suggests that there is little risk that ASF spreads beyond the fence in the east and the north. Comparing figures @ref(fig:predicted-IWB-day140-scenario2) (hunting alone) and @ref(fig:predicted-IWB-day140-scenario3) (fence alone) suggests that hunting was much more effective than the fence. Without the increased hunting there is a very likely risk of continued spread beyond the fence in the north and east -- this is most likely due to the fence not being 100% effective, and more importantly, a rapid decrease in IWB brought about by the increased hunting pressure. The posterior distribution of fence efficacy parameter \omega_\text{Fence} is summarised in table @ref(tab:posteriorOmegaFence). Clearly, there are large uncertainties concerning exact value of the fence efficacy. However, as indicated in our previous report, unless the fence is 100% efficient, an infectious wild boar can be expected to eventually carry ASF beyond the fence.

mcmcOut <- read.table(mcmcFileModel10, header=TRUE)
mySummary <- mean95(ilogit(mcmcOut$logit_omegaFence))
mySummary <- matrix(mySummary, nrow=1)
colnames(mySummary) <- c("2.5%","Mean","97.5%")
rownames(mySummary) <- ""
kable(mySummary,
      row.names = FALSE, ##  "p_~attractI~",
      format = "latex",
      caption = "Summary of the posterior distribution of fence efficacy.",
      digits = 3,
      booktabs = TRUE
      )

cap <- "Expected densities of infectious wild boar on day 140 under control scenario 2: no fence (purple dashed), increased hunting zone (brown), 1km active search (not shown)."
print(map_expectedI_scenario2_140)

cap <- "Expected densities of infectious wild boar on day 140 under control scenario 3: fence (purple) , no increased hunting zone (brown dashed), 1km active search (not shown)."
print(map_expectedI_scenario3_140)

Comparing figures @ref(fig:predicted-IWB-day140-scenario3) (fence alone) and @ref(fig:predicted-IWB-day140-scenario4) (no fence, no increased hunting) suggests that under both scenarios ASF is free to escape to the north and east. The fence slows this escape slightly (as indicated by lower IWB densities), but the effect is small and most likely insignificant in the long term.

cap <- "Expected densities of infectious wild boar on day 140 under control scenario 4: no fence (dashed purple), no increase in hunting zone (drown dashed), 1km active search (not shown)."
print(map_expectedI_scenario4_140)

Interestingly, in all four scenarios, there is an escape of ASF to the south of the hunting zone. The magnitude of this escape is much lower in the pressence of increased hunting pressure (figures @ref(fig:predicted-IWB-day140-scenario1) and @ref(fig:predicted-IWB-day140-scenario2)) however, left unchecked this new focii will be free to grow to the extent seen in the no-control scenario (figure @ref(fig:predicted-IWB-day140-scenario4))

The expected temporal dynamics of IWB density under all eight control scenarios are pressented in figure @ref(fig:IWB-scenarios). The results clearly indicate that the trajectories of these eight scenarios cluster into three groups, and that within each group the dynamics are indistinguishably similar. The most effective scenarios are 1, 2, 5 and 6, where the expected density of IWB on day 140 is r tail(IWB_scenarios,1)[1], r tail(IWB_scenarios,1)[2], r tail(IWB_scenarios,1)[5] and r tail(IWB_scenarios,1)[6] respectively. For scenarios 3 and 7 the expected density of IWB on day 140 is r tail(IWB_scenarios,1)[3] and r tail(IWB_scenarios,1)[7] respectively. And for scenarios 4 and 8 the expected density of IWB on day 140 is r tail(IWB_scenarios,1)[4] and r tail(IWB_scenarios,1)[8] respectively.

These results clearly indicate that, in our model at least:

1. increasing the radius of the active search will have close-to-zero impact;
2. fencing perhaps has some effect on slowing spread, but there is currently very little evidence suggesting that fence efficacy is high;
3. by far the most effective intervetion, that we can detect in the data, has been the increased hunting pressure;
4. without the increased hunting pressure, the epidemic would have spread beyond the fence in the north and east and the extent of spread in the south would have been far greater.
5. none of the control scenarios completely cut spread to the south of the increased hunting zone.

cap <- "Expected total infectious wild boar (IWB) by scenario. Scenarios that included increased hunting pressure (1, 2, 5 and 6) were the most effective. Scenarios with no control measures (4 and 8) were the least effective. Note, we have not presented variation about the expected value, but can do so if requested to."
nScenarios <- nrow(scenariosWBR2)
days <- as.integer(sub("D","",row.names(IWB_scenarios)))
plot(days, IWB_scenarios[,1], typ="n", xlab="Day", ylab="Total IWB", ylim=c(0,1.2*max(IWB_scenarios)))
myCols = rainbow(nScenarios)
for (ii in 1:nScenarios)
  lines(days, IWB_scenarios[,ii], col=myCols[ii])
abline(v=80,col="lightgrey")
abline(h=0)

Perhaps one reason increasing the radius of the acive search will be ineffective is because the associated detection probabilities remain rather low (Table @ref(tab:table-probAS)). Another likely reason is that living IWB are more important in the rapid transmission of ASF than infectious carcasses, as indicated by the posterior distribution of \omega_\text{AttractI} (Table @ref(tab:posteriorAttractI)), in other words, the virus spreads very rapidly among living wild boar such that the effects of infectious carcasses remaining in the landscape become negligible within the time frame of the current study.

kable(tab_probAS,
  format = "latex",
  caption = "Mean and 95 percent credibility intervals for daily active search detection probabilities.
             'Central' refers to pixels where an infected wild boar sparked an active search.
             'Neighbour' refers to adjacent pixels where the active search may also partially take place.
             1km and 2km refers to the two possible active search radii considered.",
  digits = 3,
  booktabs = TRUE
)

mcmcOut <- read.table(mcmcFileModel10, header=TRUE)
mySummary <- mean95(ilogit(mcmcOut$logit_attractI))
mySummary <- matrix(mySummary, nrow=1)
colnames(mySummary) <- c("2.5%","Mean","97.5%")
rownames(mySummary) <- ""
kable(mySummary,
      row.names = FALSE, ##  "p_~attractI~",
      format = "latex",
      caption = "Summary of the posterior distribution of the relative importance of IWB compared to carcasses in transmission events.",
      digits = 3,
      booktabs = TRUE
      )

Finally, we use our simulations to quantify the expected exposure to IWB at each of the outdoor pig farms in the affected area. We define exposure as the sum, over days 0-140, of the number of IWB within the pixel containing the pig farm. These exposures are pressented in ranked order in Table @ref(tab:table-farm-exposure-IWB-outdoor).

kable(pigSitesExposure %>% st_drop_geometry() %>%
  filter(is_outdoor==1) %>%
  select(population_id, size, FarmType, practice, exposure, tPeakExposure, knownCases),
  format = "latex",
  caption = "Table of outdoor farms most exposed to infectious wild boar, where: \\emph{exposure} is the sum (over days 0-140) of the expected number of infectious wild boar within each farm's hexagonal pixel; and tPeakExposure indicates the timing (day) of the greatest local density of infectious wild boar.",
  digits = 3,
  booktabs = TRUE
)

\clearpage

Conclusions and prediction requests

• With the current (D110) policy of disease control on farms, the disease is most likely to continue spreading locally between neighbouring farms within the fenced area and buffer zone. However, the risk of infectious animal shipments out of the area is very low according to our model (and it's assumptions).

• We identified farms with site-ids r tab_farms_pred_CSIF_D80 %>% filter(!detected, p_pred >= .3) %>% pull(id) %>% paste(collapse = ", ") as likely (p \geq .3) to have been infected in this period (D80 - D110). To a lesser extent (.1 \leq p < .3), farm sites r tab_farms_pred_CSIF_D80 %>% filter(!detected, p_pred > .1, p_pred < .3) %>% pull(id) %>% sort %>% paste(collapse = ", ") also have some risk of being infected. A few other farms have even smaller risks (see Table @ref(tab:table-farms-pred)).

• The total number of animals affected in these farms are r tab_farms_pred_CSIF_D80 %>% filter(!detected, p_pred >= .3) %>% pull(size) %>% sum in the first case and r tab_farms_pred_CSIF_D80 %>% filter(!detected, p_pred > .1, p_pred < .3) %>% pull(size) %>% sum in the second.

• The risk associated with pig trade is very low (below 1\% probability of infection for all farms).

• The expected number of affected animals, under the current control strategy (in addition to those from already detected farms) is r (exp_nanimals = exp_pigs_cs %>% filter(CS == "IF") %>% pull("Exp pigs")) %>% round() out of r round((total_animals = pig_sites_risks_D80 %>% filter(!detected) %>% pull(size) %>% sum) %>% "/"(1e6)) M heads (r round(exp_nanimals / total_animals * 100, 2)%).

• The current control strategy on farms works well. The combination of culling of infected farms and ban trading in protection and surveillance zones are really effective to contain the spread of the disease.

• Among the additional culling measures, the one with largest impact would be culling herds in the viccinity of wild boar carcasses. But the reduction is marginal (see Table @ref(tab:exp-pigs-cs)). However, these results are preliminary as we didn't fully account for the culling strategy in the assessment of the transmission via fomites. Thus, the impacts are likely a bit larger.

• The strategy consisting on increasing the size of the active-search area around wild boar carcasses from 1 to 2 km was not assessed, but we think that it would not have a major impact.

• The current strategy of maintaining an increased hunting pressure has been highly effective in reducing the size of the ASF reservoir in wild boar.

• However, we expect the disease to spread beyond the southern limit of the hunting zone before day 140. Thus additional control measures will be required to the south of the buffer zone to avoid a secondary wave after day 140.

• A number of outdoor farms have been identified as at risk from the travelling wave of ASF among wild boar. If we apply a cutoff of at exposure=5, these farms include r pigSitesExposure %>% st_drop_geometry() %>% filter(is_outdoor==1 & knownCases == FALSE & exposure > 5) %>% select(population_id) %>% pull(population_id) %>% paste(collapse = ", "). Farms with an exposure of 1-5 are listed in table @ref(tab:table-farm-exposure-IWB-outdoor). For some of these farms, the peak exposure time is estimated at 140 days (the maximum in our simulations), suggesting that a subsequent secondary wave in the south could have serious consequences if left uncontrolled.