library(tidyverse)
## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr
## Conflicts with tidy packages ----------------------------------------------
## filter(): dplyr, stats
## lag(): dplyr, stats
library(knitr)
library(beepr)
library(doParallel)
## Loading required package: foreach
##
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
##
## accumulate, when
## Loading required package: iterators
## Loading required package: parallel
library(foreach)
library(parallel)
set.seed(1994)
# Define a helper function to use in the SimEpidemic function
infect <- function(inf.status, infectivity, daily.contacts) {
infected <- which(inf.status == 1)
if(length(infected)==0) {
return(inf.status)
}
day.status <- sapply(1:length(inf.status), function(i) {
contacts <- sample.int(daily.contacts, replace = FALSE)
# Use the fact that rbinom can take a vector of probabilities
contact.results <- rbinom(daily.contacts,
1,
p = infectivity*inf.status[contacts])
return(any(contact.results==1))
})
return(ifelse(day.status, 1, 0))
}
# Define the SimEpidemic function
SimEpidemic <- function(N, init.p, infectivity, daily.contacts) {
curr.status <- rbinom(N, 1, init.p)
week.status <- NULL
for(day in 2:7) {
week.status <- rbind(week.status, curr.status)
curr.status <- infect(curr.status, infectivity, daily.contacts)
}
return(list(N.cases = rowSums(week.status),
N.infected = N - sum(colSums(week.status)==0)))
}
# Testing the code
N <- 2500
init.p <- 0.5
infectivity <- 0.1
daily.contacts <- 10
epid <- SimEpidemic(N, init.p, infectivity, daily.contacts)
#####Task 1#####
Ns <- c(100, 2500)
init.ps <- c(0.10, 0.20, 0.50)
infectivitys <- 0.10
daily.contactss <- 10
inputs <- expand.grid(Ns, init.ps, infectivitys,
daily.contactss)
### Use foreach and doParallel to do this on Windows (and on UNIX and Mac)
ncores <- detectCores()-1 #Set ncores to something less than or equal to available cores
cl <- makeCluster(ncores) # creates a cluster with <ncore> cores
registerDoParallel(cl) # register the cluster
Nsim <- 100
results <- foreach(j = 1:Nsim) %dopar% {
set.seed(j)
return(sapply(purrr::pmap(.f = SimEpidemic,
.l = list(inputs[, 1], inputs[, 2], inputs[, 3],
inputs[, 4])),function(x) x[2]))
}
stopCluster(cl) #Stop executing in parallel
Final_Results <- data.frame(N = rep(inputs[, 1], 100),
init.p = rep(inputs[, 2], 100),
N.infected = results %>% unlist())
Results <- Final_Results %>% group_by(N, init.p) %>%
summarise(Avg.N.Infected = mean(N.infected),
SD = sd(N.infected))
#####Task 3#####
Ns <- 500
init.ps <- 0.20
infectivitys <- c(0.20, 0.50, 0.80)
daily.contactss <- c(10, 25)
inputs2 <- expand.grid(Ns, init.ps, infectivitys,
daily.contactss)
Plot_Data <- pmap(.f = SimEpidemic,
.l = list(inputs2[, 1], inputs2[, 2], inputs2[, 3],
inputs2[, 4])) %>%
map(.f = function(x) x[1]) %>% bind_cols()
Plot_Data_long <- reshape(Plot_Data, idvar = "Day",
v.names = "N.cases", direction = "long",
varying = colnames(Plot_Data)) %>% rename(Scenario = time)
Plot_Data_long$Scenario <- Plot_Data_long$Scenario %>% as.factor()
new.levels <- paste0("infectivity = ", inputs2[, 3], ", \n daily contacts = ",
inputs2[, 4])
#For N = 500, init.p = 0.20
Plot_Data_long$Scenario <- Plot_Data_long$Scenario %>% recode_factor(`1` = new.levels[1], `2` = new.levels[2],
`3` = new.levels[3], `4` = new.levels[4],
`5` = new.levels[5], `6` = new.levels[6])
#Plotting the results
p <- ggplot(data = Plot_Data_long, aes(x = Day, y = N.cases, color = Scenario)) +
geom_point() + ggtitle("Number of cases by day \n N = 500, init.p = 0.20") +
guides(color = guide_legend(nrow = 3)) +
theme(legend.direction = 'horizontal',
legend.position = 'bottom',
legend.key = element_rect(size = 5),
legend.key.size = unit(1.5, 'lines'),
plot.title = element_text(hjust = 0.5)) +
geom_line()
p
Results %>% kable()
| 100 |
0.1 |
35.96 |
28.073828 |
| 100 |
0.2 |
59.75 |
24.816244 |
| 100 |
0.5 |
89.44 |
9.276145 |
| 2500 |
0.1 |
786.75 |
656.207613 |
| 2500 |
0.2 |
1450.47 |
570.047924 |
| 2500 |
0.5 |
2248.16 |
204.801886 |