diff --git a/nimble/model12-M1/simulation_model12-M1.R b/nimble/model12-M1/simulation_model12-M1.R
index 102f2ce1697b452e9be47c87ea184d2b903f91c2..48f4f49e0961a5d846c71231104334e2bf06099d 100644
--- a/nimble/model12-M1/simulation_model12-M1.R
+++ b/nimble/model12-M1/simulation_model12-M1.R
@@ -39,7 +39,7 @@ tStop = 365 ## For 30 day prediction beyond day 80
initWildboarM1 = setWildBoarModelInitialValues(Hex = hexM1,
bboxHex = bboxM1,
tStop = tStop,
- pHome = pHomeSim4km[1]/sum(pHomeSim4km)
+ pHome = pHomeSim4km
)
constWildboarM1 = setWildBoarModelConstants(Hex = hexM1,
HexCentroids = hexCentroidsM1,
@@ -122,20 +122,19 @@ nSteps1 = nSteps + 1
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
## VARIATION MEAN WITH rWBdensity
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-cWildBoarModelM1$sampleEboar<-0
-Mean<-seq(2,16,by=2)
-WB_data<-matrix(NA,nrow=nrow(hexM1),ncol=length(Mean))
-colnames(WB_data)<-paste0("Density",Mean)
-Emergence_date<-matrix(NA,nrow=nrow(hexM1),ncol=length(Mean))
-
+cWildBoarModelM1$sampleEboar <- 0
+Mean <- seq(2,16,by=2)
+WB_data <- matrix(NA,nrow=nrow(hexM1),ncol=length(Mean))
+colnames(WB_data) <- paste0("Density",Mean)
+Emergence_date <- matrix(NA,nrow=nrow(hexM1),ncol=length(Mean))
# function rWBdensity(Mean[i],SD=20,hexM1)
for (i in 1:length(Mean)){
- hexM1sim<-rWBdensity(Mean[i],SD=2,hex=hexM1)
- WB_data[,i]=hexM1sim$E_wildboar
- cWildBoarModelM1$Eboar<-hexM1sim$E_wildboar
+ hexM1sim <- rWBdensity(Mean[i],SD=2,hex=hexM1)
+ WB_data[,i] <- hexM1sim$E_wildboar
+ cWildBoarModelM1$Eboar <- hexM1sim$E_wildboar
simulate(cWildBoarModelM1, detNodes)
- simE <- t(cWildBoarModelM1$epidemic[iE, 1:nSteps1, 1:nPops])
- Emergence_date[,i]=apply(simE,1, function(x) min(which(x>0)))-1
+ simE <- t(cWildBoarModelM1$epidemic[iE, 1:nSteps1, 1:nPops])
+ Emergence_date[,i] <- apply(simE,1, function(x) min(which(x>0)))-1
}
Emergence_date[sapply(Emergence_date, is.infinite)]<- NA
colnames(Emergence_date)=paste0("Emergence_mean", Mean)
@@ -147,53 +146,65 @@ hexSim <- hexM1sim %>% select(!"E_wildboar") %>% cbind(Emergence_date) %>% cbind
#############
map1_rWBdensity<-list()
for(m in Mean){
- map1_rWBdensity[[which(m==Mean)]] <- tm_shape(hexSim) + tm_polygons(paste0("Density",m), breaks=sort(unique(c(0,1,pretty(c(1,max(WB_data)))))), legend.show = FALSE)
+ map1_rWBdensity[[which(m==Mean)]] <- tm_shape(hexSim) + tm_polygons(paste0("Density",m), breaks=sort(unique(c(0,1,pretty(c(1,max(WB_data)))))), legend.show = FALSE)
}
-map1_rWBdensity[[1+which(m==Mean)]] <- tm_shape(hexSim) + tm_polygons(paste0("Density",m),title="Density", breaks=sort(unique(c(0,1,pretty(c(1,max(WB_data))))))) +
- tm_shape(hexSim) + tm_polygons(col="white", border.col="white") + tm_layout(legend.position = c ("LEFT", "TOP"))
+map1_rWBdensity[[1+which(m==Mean)]] <- tm_shape(hexSim) + tm_polygons(paste0("Density",m),title="Density", breaks=sort(unique(c(0,1,pretty(c(1,max(WB_data))))))) +
+ tm_shape(hexSim) + tm_polygons(col="white", border.col="white") + tm_layout(legend.position = c ("LEFT", "TOP"))
-X11() # show the map in a new window
+# X11() # show the map in a new window
tmap_arrange(map1_rWBdensity)
+
+#############
+## R0 maps ##
+#############
+debug(get_R0)
+### doesn't find hex$E_wildboar within get_R0
+hexSim$R0 <- get_R0(hexSim, mu = 1/(5*365), model=cWildBoarModelM1, densityCol="Density2")
+print(tm_shape(hexSim)+tm_polygons(paste0("R0"), breaks=sort(unique(c(0,1,pretty(hexSim$R0))))))
+
+hexSim$R0_Neigh <- get_R0_Neigh(hexSim, mu = 1/(5*365),neighVec=neighVec, model=cWildBoarModelM1)
+print(tm_shape(hexSim)+tm_polygons(paste0("R0_Neigh")))
+
+
###############
# EMERGENCE MAP
###############
map2_rWBdensity<-list()
for(m in Mean){
- map2_rWBdensity[[which(m==Mean)]]<-tm_shape(hexSim) + tm_polygons(paste0("Emergence_mean",m), legend.show = FALSE)
+ map2_rWBdensity[[which(m==Mean)]]<-tm_shape(hexSim) + tm_polygons(paste0("Emergence_mean",m), legend.show = FALSE)
}
map2_rWBdensity[[1+which(m==Mean)]]<-tm_shape(hexSim) + tm_polygons(paste0("Emergence_mean",m), title = "Emergence_mean") +
tm_shape(hexSim) + tm_polygons(col="white", border.col="white") + tm_layout(legend.position = c ("LEFT", "TOP"))
-X11() # show the map in a new windo
+# X11() # Open a new window for plotting figures
tmap_arrange(map2_rWBdensity)
##%%%%%%%%%%%%%%%%%%%%
## VARIATION SD ONLY
##%%%%%%%%%%%%%%%%%%%%
-cWildBoarModelM1$sampleEboar<-0
-SD<-seq(0,10,by=2)
-SD_data<-matrix(NA,nrow=nrow(hexM1),ncol=length(SD))
-Emergence_date<-matrix(NA,nrow=nrow(hexM1),ncol=length(SD))
+cWildBoarModelM1$sampleEboar <- 0
+SD <- seq(0,10,by=2)
+SD_data <- matrix(NA,nrow=nrow(hexM1),ncol=length(SD))
+Emergence_date <- matrix(NA,nrow=nrow(hexM1),ncol=length(SD))
for (i in 1:length(SD)){
- hexM1sim<-rWBdensity(SD[i],Mean=16,hexM1)
- SD_data[,i]=hexM1sim$E_wildboar
-
- cWildBoarModelM1$Eboar<-hexM1sim$E_wildboar
+ hexM1sim <- rWBdensity(SD[i],Mean=16,hexM1)
+ SD_data[,i] <- hexM1sim$E_wildboar
+ cWildBoarModelM1$Eboar <- hexM1sim$E_wildboar
simulate(cWildBoarModelM1, detNodes)
- simE <- t(cWildBoarModelM1$epidemic[iE, 1:nSteps1, 1:nPops])
- Emergence_date[,i]=apply(simE,1, function(x) min(which(x>0)))-1
+ simE <- t(cWildBoarModelM1$epidemic[iE, 1:nSteps1, 1:nPops])
+ Emergence_date[,i] <- apply(simE,1, function(x) min(which(x>0)))-1
}
Emergence_date[sapply(Emergence_date, is.infinite)]<- NA
-colnames(Emergence_date)=paste0("Emergence_SD", SD)
-hexSim <- hexM1sim %>% select("E_wildboar") %>% cbind(Emergence_date)
+colnames(Emergence_date) <- paste0("Emergence_SD", SD)
+hexSim <- hexM1sim %>% select("E_wildboar") %>% cbind(Emergence_date)
print(tm_shape(hexSim) + tm_polygons(paste0("Emergence_SD",max(SD))))
tm_shape(hexM1sim)+tm_polygons(paste0("E_wildboar"))
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-## VARIATION MEAN WITH THE FUNCTION rWB
+## VARIATION MEAN WITH THE FUNCTION rWB
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cWildBoarModelM1$sampleEboar<-0
@@ -224,10 +235,10 @@ neighVec=attributes(hexM1)$neighVec
#############
map1_rWB<-list()
for(m in Mean){
- map1_rWB[[which(m==Mean)]] <- tm_shape(hexSim) + tm_polygons(paste0("Density",m), breaks=sort(unique(c(0,1,pretty(c(1,max(WB_data)))))), legend.show = FALSE)
+ map1_rWB[[which(m==Mean)]] <- tm_shape(hexSim) + tm_polygons(paste0("Density",m), breaks=sort(unique(c(0,1,pretty(c(1,max(WB_data)))))), legend.show = FALSE)
}
-map1_rWB[[1+which(m==Mean)]] <- tm_shape(hexSim) + tm_polygons(paste0("Density",m),title="Density", breaks=sort(unique(c(0,1,pretty(c(1,max(WB_data))))))) +
- tm_shape(hexSim) + tm_polygons(col="white", border.col="white") + tm_layout(legend.position = c ("LEFT", "TOP"))
+map1_rWB[[1+which(m==Mean)]] <- tm_shape(hexSim) + tm_polygons(paste0("Density",m),title="Density", breaks=sort(unique(c(0,1,pretty(c(1,max(WB_data))))))) +
+ tm_shape(hexSim) + tm_polygons(col="white", border.col="white") + tm_layout(legend.position = c ("LEFT", "TOP"))
X11() # show the map in a new window
tmap_arrange(map1_rWB)
@@ -237,7 +248,7 @@ tmap_arrange(map1_rWB)
###############
map2_rWB<-list()
for(m in Mean){
- map2_rWB[[which(m==Mean)]]<-tm_shape(hexSim) + tm_polygons(paste0("Emergence_mean",m), legend.show = FALSE)
+ map2_rWB[[which(m==Mean)]]<-tm_shape(hexSim) + tm_polygons(paste0("Emergence_mean",m), legend.show = FALSE)
}
map2_rWB[[1+which(m==Mean)]]<-tm_shape(hexSim) + tm_polygons(paste0("Emergence_mean",m), title = "Emergence_mean") +
tm_shape(hexSim) + tm_polygons(col="white", border.col="white") + tm_layout(legend.position = c ("LEFT", "TOP"))
@@ -301,50 +312,9 @@ plot(WB_data,Emergence_date, xlab="Wild boar density", ylab="Day of emergence",
#mu mortality rate, days
-get_R0<-function(hex, mu, model){
- tauI = model$tauI
- tauR = model$tauR
- tauC = model$tauC
- tauRot = model$tauRot
- betaI = model$betaI
- betaC = model$betaC
- pHome = model$pHome
- pAway = 1-pHome
- pHome_x_pHome = pHome*pHome
- pAway_x_pAway = pAway*pAway
- pContactIi = pHome_x_pHome + pAway_x_pAway/hex$nNeigh
- pContactCi = pHome
- R0_data <- round(hex$E_wildboar) * ( (betaI* pContactIi * tauI)/ ((tauR + tauC + mu)*(tauI + mu)) + (betaC * pContactCi * tauC * tauI) / (tauRot * (tauR + tauC + mu) * (tauI + mu)) )
- return(R0_data)
-}
hexSim$R0 <- get_R0(hexSim, mu = 1/(5*365), model=cWildBoarModelM1)
print(tm_shape(hexSim)+tm_polygons(paste0("R0"), breaks=sort(unique(c(0,1,pretty(hexSim$R0))))))
-
-get_R0_Neigh<-function(hex, mu,neighVec, model){
- tauI = model$tauI
- tauR = model$tauR
- tauC = model$tauC
- tauRot = model$tauRot
- betaI = model$betaI
- betaC = model$betaC
- pHome = model$pHome
- nNeigh=6
- pAway = 1-pHome
- pHome_x_pAway = pHome*pAway
- pAway_x_pAway = pAway*pAway
- pContactIj = pHome_x_pAway*2/nNeigh + pAway_x_pAway*2 / (nNeigh*nNeigh)
- pContactCj = pAway/nNeigh
- neighStart = hex$neighStart # Indicates start of sub-vector of neighVec for each hexagonal patch
- neighEnd = hex$neighEnd # Indicates end of sub-vector of neighVec for each hexagonal patch
- #neighVec = attributes(hex)$neighVec
- EWB_neigh <- vector(mode="numeric", length=nrow(hex))
- for (i in 1:nrow(hex)){
- EWB_neigh[i] <- sum(round(hex$E_wildboar[neighVec[neighStart[i]:neighEnd[i]]]))
- }
- R0_Neigh = EWB_neigh * ((betaI* pContactIj * tauI)/ ((tauR + tauC + mu)*(tauI + mu)) + (betaC * pContactCj * tauC * tauI) / (tauRot * (tauR + tauC + mu) * (tauI + mu)) )
- return(R0_Neigh)
-}
hexSim$R0_Neigh <- get_R0_Neigh(hexSim, mu = 1/(5*365),neighVec=neighVec, model=cWildBoarModelM1)
print(tm_shape(hexSim)+tm_polygons(paste0("R0_Neigh")))
diff --git a/src/functions.R b/src/functions.R
index a81c88931b69fcdb64f51227c9db558d918c9fef..73d02b42a870fe26e66494fb4dd56c3edd316444 100644
--- a/src/functions.R
+++ b/src/functions.R
@@ -2514,75 +2514,74 @@ plot_region_at_risk <- function(
# wild boar density log normal distribution
-rWBdensity <- function(Mean, SD, hex){
+rWBdensity <- function(Mean, SD, hex, Round=TRUE){
if (SD>0){
- Var=SD^2
+ Var = SD^2
sigma = sqrt(log(exp(log(Var)-2*log(Mean))+1))
- mu= log(Mean) - (sigma^2/2)
- hex$E_wildboar<-rlnorm(nrow(hex),meanlog=mu, sdlog=sigma )
+ mu = log(Mean) - (sigma^2/2)
+ hex$E_wildboar = rlnorm(nrow(hex),meanlog =mu, sdlog=sigma)
} else{
hex$E_wildboar = Mean
}
+ if (Round)
+ hex$E_wildboar = round(hex$E_wildboar)
hex$log10_ewb1 = log10(hex$E_wildboar +1)
return(hex)
}
-rWBdensity2 <- function(Mean, SD, hex, i){
+rWBdensity2 <- function(Mean, SD, hex, i, Round=TRUE){
if (SD>0){
Var=SD^2
sigma = sqrt(log(exp(log(Var)-2*log(Mean))+1))
mu= log(Mean) - (sigma^2/2)
- hex$E_wildboar[i]<-rlnorm(1,meanlog=mu, sdlog=sigma )
+ hex$E_wildboar[i] = rlnorm(1,meanlog=mu, sdlog=sigma )
} else{
hex$E_wildboar[i] = Mean
}
+ if (Round)
+ hex$E_wildboar[i] = round(hex$E_wildboar[i])
hex$log10_ewb1[i] = log10(hex$E_wildboar[i] +1)
return(hex)
}
rWB <- function(Mean, SD, hex, niter=0, hexCentroids, threshold=5/100){
- Var=SD^2
- sigma = sqrt(log(exp(log(Var)-2*log(Mean))+1))
- mu= log(Mean) - (sigma^2/2)
- #num<-0
- center<-identify_center(hex=hex,hexCentroids=hexCentroids)
+ Var = SD^2
+ sigma = sqrt(log(exp(log(Var)-2*log(Mean))+1))
+ mu = log(Mean) - (sigma^2/2)
+ center = identify_center(hex=hex,hexCentroids=hexCentroids)
for (i in 1:nrow(hex)) {
if (runif(1) < threshold | i==center){
hex$E_wildboar[i] = rlnorm(1,meanlog=mu, sdlog=sigma )
hex$log10_ewb1 = log10(hex$E_wildboar +1)
- #num<-num+1
}else{
- hex$E_wildboar[i] <- 0
+ hex$E_wildboar[i] = 0
}
}
- #print(num)
####
if (niter > 0){
for (iter in 1:niter){
- hexnew <- hex
+ hexnew = hex
for (i in 1:nrow(hex)) {
if (hex$E_wildboar[i]>0) {
- #num2<-num2+1
- voisin <- attributes(hex)$neighVec[hex$neighStart[i]:hex$neighEnd[i]]
+ voisin = attributes(hex)$neighVec[hex$neighStart[i]:hex$neighEnd[i]]
for (j in voisin) {
if (j >= 1 & j <= nrow(hex)) {
- hexnew$E_wildboar[j] <- hex$E_wildboar[i]
+ hexnew$E_wildboar[j] = hex$E_wildboar[i]
}
}
}
}
- hex=hexnew
- #print(tm_shape(hexnew) + tm_polygons("E_wildboar"))
+ hex = hexnew
}
}
- hex$E_wildboar=as.integer(hex$E_wildboar)
+ hex$E_wildboar = as.integer(round(hex$E_wildboar))
return(hex)
}
identify_center<-function(hex,hexCentroids){
- nPops = nrow(hex)
- xIntro=mean(st_bbox(hex)[c(1,3)])
- yIntro=mean(st_bbox(hex)[c(2,4)])
+ nPops = nrow(hex)
+ xIntro = mean(st_bbox(hex)[c(1,3)])
+ yIntro = mean(st_bbox(hex)[c(2,4)])
xCentroid = (hexCentroids %>% st_coordinates())[,"X"]
yCentroid = (hexCentroids %>% st_coordinates())[,"Y"]
distIntro = nimNumeric(length=nPops)
@@ -2592,3 +2591,50 @@ identify_center<-function(hex,hexCentroids){
hexIntro = which(distIntro[1:nPops]==min(distIntro[1:nPops]))
return(hexIntro)
}
+
+#########################
+## Mahe's R0 functions ##
+#########################
+get_R0 <- function(hex, mu, model, densityCol){
+ tauI = model$tauI
+ tauR = model$tauR
+ tauC = model$tauC
+ tauRot = model$tauRot
+ betaI = model$betaI
+ betaC = model$betaC
+ pHome = model$pHome
+ pAway = 1-pHome
+ pHome_x_pHome = pHome*pHome
+ pAway_x_pAway = pAway*pAway
+ pContactIi = pHome_x_pHome + pAway_x_pAway/hex$nNeigh
+ pContactCi = pHome
+ WBdensity = eval(parse(text=paste0("hex$",densityCol)))
+ R0_data = WBdensity * ((betaI* pContactIi * tauI)/ ((tauR + tauC + mu)*(tauI + mu)) + (betaC * pContactCi * tauC * tauI) / (tauRot * (tauR + tauC + mu) * (tauI + mu)) )
+ return(R0_data)
+}
+
+
+
+get_R0_Neigh<-function(hex, mu,neighVec, model, densityCol){
+ tauI = model$tauI
+ tauR = model$tauR
+ tauC = model$tauC
+ tauRot = model$tauRot
+ betaI = model$betaI
+ betaC = model$betaC
+ pHome = model$pHome
+ nNeigh = 6
+ pAway = 1-pHome
+ pHome_x_pAway = pHome*pAway
+ pAway_x_pAway = pAway*pAway
+ pContactIj = pHome_x_pAway*2/nNeigh + pAway_x_pAway*2 / (nNeigh*nNeigh)
+ pContactCj = pAway/nNeigh
+ neighStart = hex$neighStart # Indicates start of sub-vector of neighVec for each hexagonal patch
+ neighEnd = hex$neighEnd # Indicates end of sub-vector of neighVec for each hexagonal patch
+ EWB_neigh = vector(mode="numeric", length=nrow(hex))
+ for (i in 1:nrow(hex)){
+ EWB_neigh[i] =- sum(round(hex$E_wildboar[neighVec[neighStart[i]:neighEnd[i]]]))
+ }
+ R0_Neigh = EWB_neigh * ((betaI* pContactIj * tauI)/ ((tauR + tauC + mu)*(tauI + mu)) + (betaC * pContactCj * tauC * tauI) / (tauRot * (tauR + tauC + mu) * (tauI + mu)) )
+ return(R0_Neigh)
+}