https://github.com/morgankain/RRV_HostVectorCompetence
Tip revision: be7e87c3c4c8af0420a8dd42cdcff5586fdbad90 authored by Morgan Kain on 25 May 2021, 16:23:11 UTC
Merge pull request #1 from morgankain/add-license-1
Merge pull request #1 from morgankain/add-license-1
Tip revision: be7e87c
13_data.community.R
################################################################################
### Set up the mosquito and host communities for competence and R0 estimates ###
################################################################################
## Note: Because we do not have mosquito physiological responses and blood meals for all of the species in the community
## we simply scale up the mosquitoes that we do have these data for to assume that they are the whole mosquito community.
## For Brisbane the mosquitoes we do have data for are ~ 90% of the community, and for Cairns a bit over 70%
## Convert to proportions
mosq_prop <- mosq_prop %>% mutate(prop = count / sum(count))
## Pull out only the mosquito species that we have BOTH physiological data and mosquito biting data
mosq_prop_for_R0 <- mosq_prop %>%
filter(mos_species %in% as.character(unique(h_to_m_trans.pred$species))) %>%
filter(mos_species %in% as.character(mosq_blood.prop$mos_species))
## Re-scale
mosq_prop_for_R0$prop <- mosq_prop_for_R0$prop / sum(mosq_prop_for_R0$prop)
## All the hosts that we have titer data for
physiology_species.h <- as.character(unique(tot_titer$host))
physiology_species.m <- as.character(unique(h_to_m_emp$species))
## All the hosts that show up in the community of interest
ecological_species.h <- as.character(unique(host_prop$host_species))
ecological_species.m <- as.character(unique(mosq_blood.prop$mos_species))
## Connections here between hosts with physiological responses and ecological data are the following:
## Note, some simplifications are needed here to align the various data sources that we have
## A potentially minor problem with the simple average taken here is that titer profiles aren't getting weighted
## appropriately by abundance (e.g. within marsupials); however we lack that level of detailed abundance data
## bushtail possum -> possum
## agile wallaby -> within marsupial
## grey kangaroo -> within marsupial
## bandicoot -> dropped
## marsupial mouse -> dropped
## horse -> horse
## sheep -> sheep
## pig -> dropped because they are not allowed in urban settings
## cow -> cattle
## grey h f fox -> flying fox
## rat -> rat
## rabbit -> rabbit
## chicken -> within bird
## black duck -> within bird
## little corella -> within bird
## human -> human
#####
## Steps to adjust h_to_m_trans, which already has titer incorporated
#####
## 1) Select only the mosquitoes in the community for which we have biting preferences
mosq.which <- which(dimnames(h_to_m_trans_all_samps)[[3]] %in% ecological_species.m)
mosq.ordered <- dimnames(h_to_m_trans_all_samps)[[3]][mosq.which]
mosq_prop_for_R0 <- mosq_prop_for_R0[match(mosq.ordered, as.character(mosq_prop_for_R0$mos_species)), ]
## 2) Average marsupials (actually macropods given the removal of bandicoot and marsupial mouse)
marsupials <- c("agile_wallaby", "grey_kangaroo")
marsupials.data <- h_to_m_trans_all_samps[, which(dimnames(h_to_m_trans_all_samps)[[2]] %in% marsupials), , ]
marsupials.data.s <- apply(marsupials.data, c(1, 3, 4), mean)
marsupials.data.s <- marsupials.data.s[, mosq.which, ]
## Do the same manipulation for AUC for the first panel of Figure 2: host. Not used in the rest of the model
marsupials.data.AUC <- host_titer_AUC_all_samps[which(dimnames(host_titer_AUC_all_samps)[[1]] %in% marsupials), ]
marsupials.data.s.AUC <- colMeans(marsupials.data.AUC)
## 3) Average birds
birds <- c("chicken", "black_duck", "little_corella")
birds.data <- h_to_m_trans_all_samps[, which(dimnames(h_to_m_trans_all_samps)[[2]] %in% birds), , ]
birds.data.s <- apply(birds.data, c(1, 3, 4), mean)
birds.data.s <- birds.data.s[, mosq.which, ]
## Do the same manipulation for AUC for the first panel of Figure 2: host. Not used in the rest of the model
birds.data.AUC <- host_titer_AUC_all_samps[which(dimnames(host_titer_AUC_all_samps)[[1]] %in% birds), ]
birds.data.s.AUC <- colMeans(birds.data.AUC)
h_to_m_trans_adj_to_com <- array(dim = c(length(inf_days)
## total length - summaries + cat and dog which don't have data because they never show titer
, length(ecological_species.h)
, length(mosq.which)))
h_to_m_trans_all_samps_adj_to_com <- array(dim = c(length(inf_days)
## total length - summaries + cat and dog which don't have data because they never show titer
, length(ecological_species.h)
, length(mosq.which), n_samps))
dimnames(h_to_m_trans_adj_to_com)[[2]] <- ecological_species.h
dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] <- ecological_species.h
## Create a smaller matrix for AUC as well
host_titer_AUC_all_samps_adj_to_com <- matrix(nrow = length(ecological_species.h), ncol = n_samps, data = 0)
dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] <- ecological_species.h
### And then subset all of the mosquito data to these species
dimnames(h_to_m_trans_all_samps_adj_to_com)[[3]] <- mosq.ordered
## The summaries and cat and dog
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "bird", , ] <- birds.data.s
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "macropod", , ] <- marsupials.data.s
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "cat", , ] <- 0
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "dog", , ] <- 0
## The rest of the species (select the estimates for all the mosquitoes we are retaining and make the name change if necessary)
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "human", , ] <-
h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "human", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "possum", , ] <-
h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "brushtail_possum", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "flying_fox", , ] <-
h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "grey_h_f_fox", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "cattle", , ] <-
h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "cow", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "horse", , ] <-
h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "horse", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "sheep", , ] <-
h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "sheep", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "rabbit", , ] <-
h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "rabbit", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "rat", , ] <-
h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "rat", mosq.which, ]
## fill in for AUC
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "bird", ] <- birds.data.s.AUC
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "macropod", ] <- marsupials.data.s.AUC
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "cat", ] <- 0
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "dog", ] <- 0
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "human", ] <-
host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "human", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "possum", ] <-
host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "brushtail_possum", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "flying_fox", ] <-
host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "grey_h_f_fox", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "cattle", ] <-
host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "cow", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "horse", ] <-
host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "horse", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "sheep", ] <-
host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "sheep", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "rabbit", ] <-
host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "rabbit", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "rat", ] <-
host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "rat", ]
#####
## And subset prop_inf for the community of interest
#####
prop_inf_for_R0 <- prop_inf
prop_inf_for_R0 <- prop_inf_for_R0 %>% mutate(
host = mapvalues(host
, from = c(
"brushtail_possum"
, "agile_wallaby"
, "grey_kangaroo"
, "bandicoot"
, "marsupial_mouse"
, "horse"
, "sheep"
, "pig"
, "cow"
, "grey_h_f_fox"
, "rat"
, "rabbit"
, "chicken"
, "black_duck"
, "little_corella"
, "human"
, "cat"
, "dog"
) , to = c(
"possum"
, "macropod"
, "macropod"
, "marsupial"
, "marsupial_mouse"
, "horse"
, "sheep"
, "pig"
, "cattle"
, "flying_fox"
, "rat"
, "rabbit"
, "bird"
, "bird"
, "bird"
, "human"
, "cat"
, "dog"
)
))
## Last adjustment for proportion that become infected from the infection experiments
prop_inf_for_R0 <- prop_inf_for_R0 %>%
group_by(host) %>%
summarize(tot_inf = sum(tot_inf), tot_vir = sum(tot_vir)) %>%
mutate(num_inf = tot_vir / tot_inf)
prop_inf_for_R0[prop_inf_for_R0$host == "cat" | prop_inf_for_R0$host == "dog", 4] <- 0
prop_inf_for_R0 <- prop_inf_for_R0[match(host_prop$host_species, prop_inf_for_R0$host), ]
## Grab the info for host seroprevalence for Brisbane for use in both Brisbane and Cairns since we do not have this data for Cairns
host_sero <- host_sero[match(prop_inf_for_R0$host, host_sero$host), ]
#################################################################
### Supplemental Titer plot #2: averaged titer curves (to be made, maybe...??)
#####
####
## Plot the probabilities of infection over time from the titre curves and the mosquito infection probabilities scaled
## to consider the proportion of infected hosts that show a viraemic response
####
h_to_m_trans_all_samps_adj_to_com.gg <- melt(h_to_m_trans_all_samps_adj_to_com)
names(h_to_m_trans_all_samps_adj_to_com.gg) <- c("day", "host", "mosquito", "sample", "probability")
h_to_m_trans_all_samps_adj_to_com.gg <- h_to_m_trans_all_samps_adj_to_com.gg %>%
group_by(host, day) %>%
summarise(
lwr = quantile(probability, 0.20)
, est = quantile(probability, 0.500)
, upr = quantile(probability, 0.80)
)
inf_prob_plots <- ggplot(h_to_m_trans_all_samps_adj_to_com.gg, aes(day, est)) +
geom_line() +
geom_ribbon(aes(x = day, ymin = lwr, ymax = upr), alpha = 0.25) +
facet_wrap(~host)