From 0d4cb2ced5555ff927f68954eda0c5db8f540e29 Mon Sep 17 00:00:00 2001
From: david <david.pleydell@inrae.fr>
Date: Mon, 23 Nov 2020 00:26:14 +0100
Subject: [PATCH] Recovering nearly-lost work.
---
nimble/model9/fitWildBoarModel9.R | 2 +-
nimble/model9/model_definition.R | 8 +-
src/functions.R | 40 +++++++-
src/plan.R | 58 ++++++++++-
src/simulationModel.Rmd | 5 +-
src/submission_2.Rmd | 154 +++++++++++++++++++-----------
6 files changed, 199 insertions(+), 68 deletions(-)
diff --git a/nimble/model9/fitWildBoarModel9.R b/nimble/model9/fitWildBoarModel9.R
index 80bfb0d..7526683 100644
--- a/nimble/model9/fitWildBoarModel9.R
+++ b/nimble/model9/fitWildBoarModel9.R
@@ -1,6 +1,6 @@
# source(here::here("nimble/model9/fitWildBoarModel9.R"))
# rm(list=ls())
-iNSim = 5
+iNSim = 4
source(here::here("src/packages.R"))
source(here::here("src/functions.R"))
diff --git a/nimble/model9/model_definition.R b/nimble/model9/model_definition.R
index 1336f8f..31af076 100644
--- a/nimble/model9/model_definition.R
+++ b/nimble/model9/model_definition.R
@@ -33,14 +33,18 @@ bugsBoarCode = nimbleCode({
logit(scTIntro) ~ dLogitBeta(7, 7) # Introduction date (on 0:1)
logit(pIntroX) ~ dLogitBeta(1, 1) # Introduction location relative to bounding box of whole island
logit(pIntroY) ~ dLogitBeta(1, 1) # Introduction location relative to bounding box of whole island
- logit(pHome) ~ dLogitBeta(2, 2) # Connectivity: pHome=p(wild boar in home cell). (1-pHome)=p(wild boar in neighbouring cell)
+# logit(pHome) ~ dLogitBeta(2, 2) # Connectivity: pHome=p(wild boar in home cell). (1-pHome)=p(wild boar in neighbouring cell)
logit(attractI) ~ dLogitBeta(2, 2) # Relative contact rate, i.e. how attractive is an infectious individual relative to a carcass.
logit(omegaFence) ~ dLogitBeta(2, 2) # Fence efficacy
log(tauP) ~ dnorm(0, sd=3) # Carcass detection rate - passive detection. TESTING THIS PRIOR TO AVOID DISAPPEARING TO ZERO
log(tauA) ~ dLogExp(rate = 1/1e+07) # Carcass detection rate - active detection
log(tauHArmy) ~ dLogExp(rate = 1/1e+07) # Increase in hunting rate around fenced area after day 60
log(tauPArmy) ~ dLogExp(rate = 1/1e+07) # Increase in passive search rate around fenced area after day 60
- log(beta) ~ dLogExp(rate = 1) # The scale of transmission - multiply with the realtive attractivity of Infectious and Carcasses.
+# log(beta) ~ dLogExp(rate = 1) # The scale of transmission - multiply with the realtive attractivity of Infectious and Carcasses.
+ # Changes applied to MCMC4 15h Sunday to attempt reducing beta~pHome redundancy
+ log(beta) ~ dLogExp(rate = 0.1)
+ logit(pHome) ~ dLogitBeta(10, 10)
+ ##
betaI <- beta * attractI
betaC <- beta * (1 - attractI)
tIntro <- tMin - scTIntro * tMin
diff --git a/src/functions.R b/src/functions.R
index 83f7c17..be891ab 100644
--- a/src/functions.R
+++ b/src/functions.R
@@ -66,7 +66,7 @@ aggregateWildBoarDataFast <- function(AllBoarData, Hex, stepsPerChunk=1) {
boarData <- AllBoarData %>%
filter(DET!="NAS") %>%
select(!c(DATE.CONF, HOST)) %>% ## DATE.CULLING, ID,
- filter(!is.na(rowHex)) %>% ## DP: there are some NAs around the coast where hexagons and points do not over lap
+ filter(!is.na(rowHex)) %>% ## DP: there are some NAs around the coast where hexagons and points do not over lap
select(DET, DATE.SUSP, rowHex) %>%
st_drop_geometry()
## %>% filter(rowHex <= 3)
@@ -1541,3 +1541,41 @@ firstSampleAboveMean <- function(x) {
}
return(ii)
}
+
+
+scenariosWildBoarReport2 <- function() {
+ flagFence = c(1, 0) ## if 0 -> set logit_omegaFence to -Inf
+ flagHuntZone = c(1, 0) ## if 0 -> set log_tauHArmy to -Inf
+ radiusAS = c(1, 2) ## controls prob of AS in neighbouring cells
+ scenarios = expand.grid(flagFence=flagFence, flagHuntZone=flagHuntZone, radiusAS=radiusAS)
+ return(scenarios)
+}
+
+
+getFarmExposureToInfectiousWildBoar <- function(simIWB_scenario, pigSites, HexAll, Outbreaks, thresh=1) {
+ ## HexAll row index for each pig farm
+ (pigSites$id_HexAll <- as.numeric(st_within(pigSites, HexAll)))
+ ## HexAll row index for each farm
+ (id_HexAllPolygonsWithFarms <- sort(unique(as.numeric(st_within(pigSites, HexAll)))))
+ ## HexAll row index for hexagons with expected cumulative exposure of at least one
+ (id_HexAllCasesInWB <- which(colSums(simIWB_scenario) > thresh))
+ ## Indices of hexagons with exposed farms
+ (id_HexAllWithExposedFarms <- intersect (id_HexAllPolygonsWithFarms, id_HexAllCasesInWB))
+ ## Ranked matrix of most exposed hexagons
+ (HexAllExposed <- t(t(sort(colSums(simIWB_scenario)[id_HexAllWithExposedFarms], decreasing = TRUE))))
+ (HexAllExposed <- cbind(as.integer(sub("hex","",rownames(HexAllExposed))), HexAllExposed))
+ (colnames(HexAllExposed) <- c("idHexAll", "exposure"))
+ HexAllExposed
+ ## Subset of farms and their exposure level
+ (pigSitesAtRisk <- pigSites[is.element(pigSites$id_HexAll, HexAllExposed[,"idHexAll"]),])
+ (pigSitesAtRisk$exposure <- HexAllExposed[ match(pigSitesAtRisk$id_HexAll, HexAllExposed[,"idHexAll"]),"exposure"])
+ (pigSitesAtRisk <- pigSitesAtRisk[order(pigSitesAtRisk$exposure, decreasing=TRUE),])
+ ## identify if they already had cases
+ (id_pig_cases = Outbreaks %>% filter(!is.na(ID)) %>% select(ID) %>% st_drop_geometry())
+ ## Add flag for known cases
+ pigSitesAtRisk$knownCases <- is.element(pigSitesAtRisk$population_id, id_pig_cases$ID)
+ ## Time (day) of peak exposure
+ pigSitesAtRisk$tPeakExposure = apply(simIWB_scenario[, pigSitesAtRisk$id_HexAll], 2, function(x) which(x==max(x))-1)
+ ## Output
+ pigSitesAtRisk
+}
diff --git a/src/plan.R b/src/plan.R
index 7608163..1108579 100644
--- a/src/plan.R
+++ b/src/plan.R
@@ -730,8 +730,7 @@ plan <- drake_plan(
) %>%
st_drop_geometry() %>%
select(id = population_id, detected, p_pred = p_fomites) %>%
- ## The accuracy of the risk assessment under this threshold is
- ## ridiculous. Better neglect it.
+ ## The accuracy of the risk assessment under this threshold is ridiculous. Better neglect it.
filter(p_pred > 0.001),
.id = "pathway"
) %>%
@@ -783,7 +782,60 @@ plan <- drake_plan(
## WILD BOAR SIMULATION MODEL ####
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- ## source(here("nimble/model8/fitWildBoarModel8.R"))
+ ## here("nimble/model1/fitWildBoarModel1.R") ## Estimation for report 1
+ ## here("nimble/simWildBoarFromMCMC.R") ## Simulation for report 1
+
+ ## here("nimble/model9/fitWildBoarModel9.R") ## Estimation for report 2
+ ## here("nimble/model9/simWildBoarFromMCMC9.R") ## Simulation for report 2
+
+ simDirectoryReport2 = here("nimble/model9/MCMC3/sims") # Hopefully we can "upgrade" to MCMC4/sims beofre submission
+ ,
+
+ scenariosWBR2 = scenariosWildBoarReport2()
+ ,
+
+ simIWB_meanI_111 = read.table(
+ file_in(!!here("nimble/model9/MCMC3/sims/Mean_Fence1_Army1_Radius1_nSim5000.csv") ),
+ row.names = 1, col.names = paste0("hex",0:nrow(hexAll))
+ ## !!here("data/TimeSeries_day_80_all.csv"))
+ ),
+
+ ## ## hexAll row index for each pig farm
+ ## (pig_sites$id_hexAll <- as.numeric(st_within(pig_sites, hexAll)))
+ ## ## POINTS of outdoor pig farms
+ ## (pt_outdoorSites <- pig_sites%>%filter(is_outdoor==1))
+ ## ## hexAll row index for outdoor sites
+ ## (id_hexAllPolygonsWithOutdoorFarms <- sort(unique(as.numeric(st_within(pt_outdoorSites, hexAll)))))
+ ## ## hexAll row index for hexagons with expected cumulative exposure of at least one
+ ## (id_hexAllCasesInWB <- which(colSums(simIWB_meanI_111) > 1))
+ ## ## Indices of hexagons with exposed outdoor farms
+ ## (id_hexAllWithExposedFarms <- intersect (id_hexAllPolygonsWithOutdoorFarms, id_hexAllCasesInWB))
+ ## (hexAllExposed <- t(t(sort(colSums(simIWB_meanI_111)[id_hexAllWithExposedFarms], decreasing = TRUE))))
+ ## (hexAllExposed <- cbind(as.integer(sub("hex","",rownames(hexAllExposed))), hexAllExposed))
+ ## (colnames(hexAllExposed) <- c("idHexAll", "exposure"))
+ ## hexAllExposed
+ ## (pigSitesExposure <- pig_sites[is.element(pig_sites$id_hexAll, hexAllExposed[,"idHexAll"]),])
+ ## (pigSitesExposure$exposure <- hexAllExposed[ match(pigSitesExposure$id_hexAll, hexAllExposed[,"idHexAll"]),"exposure"])
+ ## (pigSitesExposure <- pigSitesExposure[order(pigSitesExposure$exposure, decreasing=TRUE),])
+ ## (pigSitesExposure %>% st_drop_geometry)
+ ## ## identify if they already had cases
+ ## (id_pig_cases = outbreaks %>% filter(!is.na(ID)) %>% select(ID) %>% st_drop_geometry())
+ ## ## Add flag for known cases
+ ## pigSitesExposure$knownCases <- is.element(pigSitesExposure$population_id, id_pig_cases$ID)
+ ## ## Time (day) of peak exposure
+ ## pigSitesExposure$tPeakExposure = apply(simIWB_meanI_111[, pigSitesExposure$id_hexAll], 2, function(x) which(x==max(x))-1)
+
+ pigSitesExposure = getFarmExposureToInfectiousWildBoar(simIWB_meanI_111, pig_sites, hexAll, outbreaks, thresh=1)
+ ,
+
+ ## pigSitesExposure %>% filter(is_outdoor==1) %>% select(id_hexAll) %>% st_drop_geometry() %>% unlist() -> idExposedHex
+ ## x <- as.integer(sub("D","",rownames(simIWB_meanI_111)))
+ ## plot(x, simIWB_meanI_111[,1], ylim=c(0,max(simIWB_meanI_111)), typ="n", ylab="E[Infectious Wild Boar]", xlab="Day")
+ ## for (hx in idExposedHex) { # 1:nrow(hexAll)
+ ## lines(x, simIWB_meanI_111[,hx])
+ ## #which(simIWB_meanI_111[,hx]==(simIWB_meanI_111[,hx]))
+ ## }
+
##%%%%%%%%%%%%%%%%%%
## OTHER MODELS ####
diff --git a/src/simulationModel.Rmd b/src/simulationModel.Rmd
index 454a92a..03bfcc2 100644
--- a/src/simulationModel.Rmd
+++ b/src/simulationModel.Rmd
@@ -327,10 +327,9 @@ In practice, the weighted average is dominated by a very small number of simulat
Thus the weighted average looks very similar to the most likely simulation -- the main difference being that in the weighted average we see a greater number of pixels
with some (most likely very small) value between zero and one.
-[I can't get drake to embed the video (hosted on Dropbox), but you can see it by clicking here.](https://www.dropbox.com/home/challenge?preview=combined1_iRow5866_-689115.mp4)
-
-[And here is a cmplimentary video showing the weighted standard deviation.](https://www.dropbox.com/s/oboo6f8eur0jxan/standDevI_iRow5866_-689115.mp4?dl=0)
+[I can't get drake to embed the video (hosted on Dropbox), but you can see it by clicking here.](https://www.dropbox.com/home/challenge-ASF?preview=combined1_iRow5866_-689115.mp4)
+[And here is a complimentary video showing the weighted standard deviation.](https://www.dropbox.com/home/challenge-ASF?preview=standDevI_iRow5866_-689115.mp4)
diff --git a/src/submission_2.Rmd b/src/submission_2.Rmd
index 2beeba0..826cd6b 100644
--- a/src/submission_2.Rmd
+++ b/src/submission_2.Rmd
@@ -66,8 +66,10 @@ loadd(
pig_sites_pred_D50,
pig_sites_pred_D80,
pig_sites_risks_D80,
+ pigSitesExposure,
tab_farms_pred_D50,
tab_farms_pred_D80,
+ simIWB_meanI_111,
tmAdminPigTypesAndOutbreaks,
tmProportionAgricultureMap,
tmProportionForestCoverMap,
@@ -82,30 +84,30 @@ loadd(
```{r outbreaks_predD50_vs_realised}
-outbreaks_predD50_vs_realised <-
+outbreaks_predD50_vs_realised <-
full_join(
outbreaks_at_D80 %>%
- filter(HOST == "pig herd") %>%
- st_drop_geometry() %>%
+ filter(HOST == "pig herd") %>%
+ st_drop_geometry() %>%
mutate(Infection = ifelse(DATE.SUSP <= 50, "Initial", "New"))
,
- tab_farms_pred_D50 %>%
- select(-size) %>%
- mutate(predicted = p_pred > 0.05) %>%
+ tab_farms_pred_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)) %>%
+ ) %>%
+ mutate(Infection = fct_inorder(Infection)) %>%
left_join(
pig_sites,
by = c(ID = "population_id")
- ) %>%
+ ) %>%
st_sf()
```
@@ -116,11 +118,11 @@ cap <- "Predicted risks (using updated methods) and realised infections in pig s
tm_shape(
admin,
bbox = st_bbox(
- outbreaks_predD50_vs_realised %>%
+ outbreaks_predD50_vs_realised %>%
st_buffer(dist = 5e4)
)
) +
- tm_layout(bg.color = "lightblue") +
+ tm_layout(bg.color = "lightblue") +
# tm_grid(col = "grey", labels.size = 1) +
tm_scale_bar(
breaks = c(0, 50, 100, 150),
@@ -145,9 +147,9 @@ tm_shape(
# tm_dots(col = "FarmType", size=0.3, palette = "viridis") +
tm_shape(
animal_movements %>%
- filter(date >= 50 - detection_delay, date <= 50) %>%
+ filter(date >= 50 - detection_delay, date <= 50) %>%
filter(
- source %in%
+ source %in%
filter(pig_sites_pred_D50, susceptible_wb, !detected)$population_id,
dest %in% farm_pred_inf_ship_D50$dest
)
@@ -161,13 +163,13 @@ tm_shape(
```
```{r unpredicted-sites}
-unpredicted_infections <-
- outbreaks_predD50_vs_realised %>%
+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) %>%
+unpredicted_infections %>%
+ st_drop_geometry() %>%
+ select(ID, DET, starts_with("DATE"), FarmType, practice, size) %>%
kable(
format = "latex",
booktabs = TRUE,
@@ -179,7 +181,7 @@ unpredicted_infections %>%
```{r susp-contact-case}
-susp_contact_case <-
+susp_contact_case <-
filter(pig_sites_pred_D50, susceptible_wb)[
st_nearest_feature(
unpredicted_infections,
@@ -187,18 +189,18 @@ susp_contact_case <-
),
]
-dist_to_susp_contact_case <-
+dist_to_susp_contact_case <-
st_distance(
unpredicted_infections,
susp_contact_case
) %>% as.numeric
-pig_sites_pred_D80 %>%
+pig_sites_pred_D80 %>%
inner_join(
susp_contact_case %>% st_drop_geometry() %>% select(population_id)
- ) %>%
- st_drop_geometry() %>%
- select(Id = population_id, FarmType, practice, size, "day symptoms" = date_susp) %>%
+ ) %>%
+ st_drop_geometry() %>%
+ select(Id = population_id, FarmType, practice, size, "day symptoms" = date_susp) %>%
kable(
format = "latex",
booktabs = TRUE,
@@ -212,7 +214,7 @@ We did not identify it because it was an __in-door__ site which we had excluded
previous analysis.
However, it is in close proximity (~`r round(dist_to_susp_contact_case, -1)` m) to another infected very small __outdoor__ finishing farm (Table \@ref(tab:susp-contact-case))
-This suggests a fomite transmission pathway, related to proximity to infected
+This suggests a fomite transmission pathway, related to proximity to infected
farms. We include this pathway in the current prediction iteration.
\clearpage
@@ -229,7 +231,7 @@ We consider the risk higher when infected farms are in close proximity and thus
Since infected farms in the viccinity could have not yet been detected, we consider all neighbouring farms weighted by their exposition to infectious wild boar (TODO: probability of infection due to WB?)
-Another change from the previous submission is a correction of the estimation of the probability of infection via animal shipment.
+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 the destination farm site would have been detected by now.
@@ -259,20 +261,20 @@ sites_area <- pig_sites_pred_D80 %>% filter(expo_inf_wb > 0) # not only outdoo
```{r}
-distance_to_nearest_infected_site <-
+distance_to_nearest_infected_site <-
st_distance(
- pig_sites_risks_D80 %>%
+ pig_sites_risks_D80 %>%
filter(!detected, susceptible_wb),
- pig_sites_risks_D80 %>%
+ pig_sites_risks_D80 %>%
filter(detected == 1, susceptible_wb),
- ) %>%
+ ) %>%
apply(1, min)
```
```{r risk-fomites, include = FALSE, eval = FALSE}
## Assume a exponentially decreasing exposition kernel
-## with mean of 700 m (the only suspected case) and multiply
+## with mean of 700 m (the only suspected case) and multiply
## by the probability of wild boar infection, taking 1
## for already infected farms.
@@ -282,8 +284,8 @@ distance_to_nearest_infected_site <-
pig_sites_pred_D80 %>%
mutate(
p_fomites = prob_fomites(p_inf, unclass(st_distance(.)))
- ) %>%
- filter(p_fomites > 0.005) %>%
+ ) %>%
+ filter(p_fomites > 0.005) %>%
ggplot(aes(p_inf, p_fomites)) +
geom_point() +
geom_point(
@@ -299,11 +301,11 @@ pig_sites_pred_D80 %>%
# scale_y_log10() +
coord_fixed()
-idx <-
+idx <-
pig_sites_pred_D80 %>%
mutate(
p_fomites = prob_fomites(p_inf, unclass(st_distance(.)))
- ) %>%
+ ) %>%
slice_max(p_fomites, n = 5)
```
@@ -324,11 +326,11 @@ 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 <-
+infected_sites_region <-
outbreaks_at_D80 %>%
- filter(HOST == "pig herd") %>%
- select(ID) %>%
- st_set_agr("identity") %>%
+ filter(HOST == "pig herd") %>%
+ select(ID) %>%
+ st_set_agr("identity") %>%
st_intersection(st_union(wb_case_density_D80))
ggplot() +
@@ -396,21 +398,21 @@ ggplot() +
cap <- "Probability of farm infection and detection (curve), estimated from pressence/abscence data (points),
for outdoor farms in the region at risk, as a function of the exposition level."
pig_sites_risks_D80 %>%
- filter(susceptible_wb) %>%
+ filter(susceptible_wb) %>%
mutate(
detected = as.numeric(detected)
) %>%
left_join(
outbreaks_at_D80 %>%
- filter(HOST == "pig herd") %>%
- st_drop_geometry() %>%
+ 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")) %>%
+ ) %>%
+ replace_na(list(Infection = "Not infected")) %>%
ggplot(aes(expo_inf_wb, detected)) +
geom_point(aes(group = population_id, colour = Infection)) +
geom_line(
@@ -458,8 +460,8 @@ plot_area <- st_bbox(
c(
pig_sites_pred_D80 %>% filter(susceptible_wb) %>% st_geometry(),
pig_sites_pred_D80 %>% filter(susceptible_wb) %>% st_geometry() %>% `+`(c(5e4, 0))
- ) %>%
- st_union() %>%
+ ) %>%
+ st_union() %>%
st_buffer(dist = 1e4)
)
@@ -513,7 +515,7 @@ tmAdmin +
tm_shape(
animal_movements %>%
filter(
- source %in%
+ source %in%
filter(pig_sites_pred_D80, susceptible_wb)$population_id,
dest %in% farm_pred_inf_ship_D80$dest)
) +
@@ -531,7 +533,7 @@ tmAdmin +
```{r prob-infection-fomites, fig.cap = cap}
cap <- "Probability of detection due to exposition to fomites during the prediction period (D65 - D95)."
-sites_fomites <-
+sites_fomites <-
pig_sites_pred_D80 %>%
mutate(
p_fomites = prob_fomites(p_inf, unclass(st_distance(.)))
@@ -542,8 +544,8 @@ plot_area <- st_bbox(
c(
sites_fomites %>% st_geometry(),
sites_fomites %>% st_geometry() %>% `+`(c(7e4, 0))
- ) %>%
- st_union() %>%
+ ) %>%
+ st_union() %>%
st_buffer(dist = 1e4)
)
@@ -577,7 +579,7 @@ tm_shape(admin, bbox = plot_area) +
```{r table-farms-pred}
kable(
tab_farms_pred_D80 %>%
- filter(p_pred > 0.05, !detected) %>%
+ 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",
@@ -592,9 +594,9 @@ kable(
```{r prob-infection-all, fig.cap = cap}
cap <- "Probability of infection and detection by any possible route during the prediction period (D65 - D95)."
-sites_all <-
- pig_sites_pred_D80 %>%
- select(-detected, -size, -p_pred) %>%
+sites_all <-
+ pig_sites_pred_D80 %>%
+ select(-detected, -size, -p_pred) %>%
inner_join(
bind_rows(
tab_farms_pred_D80 %>%
@@ -609,8 +611,8 @@ plot_area <- st_bbox(
c(
sites_all %>% st_geometry(),
sites_all %>% st_geometry() %>% `+`(c(7e4, 0))
- ) %>%
- st_union() %>%
+ ) %>%
+ st_union() %>%
st_buffer(dist = 1e4)
)
@@ -649,6 +651,42 @@ tm_shape(admin, bbox = plot_area) +
# Wild boar
+
+```{r table-farm-exposure-IWB-outdoor}
+kable(
+ pigSitesExposure %>% st_drop_geometry() %>% # filter(knownCases == FALSE) %>%
+ filter(is_outdoor==1) %>%
+ select(population_id, size, FarmType, practice, exposure, tPeakExposure, knownCases) ,
+ caption = "Table of outdoor farms most exposed to infectious wild boar, where:
+ exposure_ is the sum, over days 0-140, of the expected number of infectious wild boar within each farm's hexagonal pixel;
+ and tPeakExposure provides the timing (day) with the gretest local density of infectious wild boar.",
+ format = "latex",
+ digits = 3,
+ booktabs = TRUE
+)
+```
+
+
+```{r table-farm-exposure-IWB-indoor}
+kable(
+ pigSitesExposure %>% st_drop_geometry() %>% # filter(knownCases == FALSE) %>%
+ filter(is_outdoor==0) %>%
+ select(population_id, size, FarmType, practice, exposure, tPeakExposure, knownCases)
+ ,
+ caption = "Table of indoor farms located in pixels with greatest exposed to infectious wild boar,
+ where:
+ exposure_ is the sum, over days 0-140, of the expected number of infectious wild boar within each farm's hexagonal pixel;
+ and tPeakExposure provides the timing (day) with the gretest local density of infectious wild boar.",
+ format = "latex",
+ digits = 3,
+ booktabs = TRUE
+)
+```
+
+
+
+
+
\clearpage
# Conclusions and prediction requests
@@ -661,5 +699,5 @@ infected in this period. To a lesser extent ($.1 \leq p \leq .35$), farm sites `
The total number of animals affected in these farms are
`r tab_farms_pred_D80 %>% filter(!detected, p_pred >= .35) %>% pull(size) %>% sum` in the first case and `r tab_farms_pred_D80 %>% filter(!detected, p_pred > .1, p_pred < .35) %>% pull(size) %>% sum` in the second.
-
+
The expected number of affected animals (in addition to those from already detected farms) is `r (exp_nanimals = tab_farms_pred_D80 %>% filter(!detected) %>% with(., sum(p_pred * size))) %>% 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)`%).
--
GitLab