#******************************* SERENGETI UNCERTAINTY ANALYSIS ***************************# # # # # #******************************************************************************************# #*************************************** PARAMETERS ********************************************# clock.start <- Sys.time() #start timing trials <- 10000 #set number of runs here days <- 3650 #days of each run start.up <- 3000 #demographic initialization time rv <- 1 pb <- txtProgressBar(title = "progress bar", min = 0, max = trials, width = 300, style = 3) sex.ratio <- rnorm(trials, 0.43, 0.02) #from Hampson 2009 litter.freq <- rnorm(trials, 0.84, .0255) litter.size <- rnorm(trials, 4.76, .153) a1 <- sex.ratio*litter.freq*litter.size #dog annual birth rate a.vec <- a1/365 #dog daily birth rate b1 <- rnorm(trials, 0.45, 0.02) #dog annual death rate b.vec <- (b1)/365 #dog daily death rate c.vec <- rnorm(trials, mean(a.vec), sd(a.vec)) #wildlife birth rate d.vec <- rnorm(trials, mean(b.vec), sd(b.vec)) #wildlife death rate x <- 5/3373 add <- rep(x/365, days) sigma.vec <- 1/rnorm(trials, 22.3, 1.28) #1/incubation period alpha.vec <- 1/rnorm(trials, 3.1, .128) #1/infectious period Kd <- 10 #carrying capacity of dogs Kw <- 3 #carrying capacity for wildlife SD.k11.vec <- rtruncnorm(trials, a=0, b=Inf, 1.01278771, 0.03767664) #dog to dog SD.k21.vec <- rtruncnorm(trials, a=0, b=Inf, 0.09145798, 0.01036436) #wildlife to dog SD.k12.vec <- rtruncnorm(trials, a=0, b=Inf, 0.95180666, 0.12316087) #dog to wildlife SD.k22.vec <- rtruncnorm(trials, a=0, b=Inf, 0.23178093, 0.05684611) nu.vec <- seq(0, 1, by = 0.1) #vaccination coverage levels record.days <- seq(1, 3650, by = 50) #time points to be saved year.marks <- seq(365, 3650, by = 365) SD.rabid <- rep(0, times = length(nu.vec)*trials*10) dim(SD.rabid) <- c(trials, length(nu.vec), 10) #array of rabid dogs by year of program SD.vaccinated <- rep(0, times = length(nu.vec)*trials*10) dim(SD.vaccinated) <- c(trials, length(nu.vec), 10) #array of vaccinated dogs by year of program #arrays for transmission throughout time #rows = trials, columns = vacc coverage, dim3 = time SD.rabid <- rep(0, times = length(nu.vec)*trials*10*5) dim(SD.rabid) <- c(trials, 5, length(nu.vec), 10) #array of rabid dogs by year of program SD.vaccinated <- rep(0, times = length(nu.vec)*trials*10*5) dim(SD.vaccinated) <- c(trials, 5, length(nu.vec), 10) #array of vaccinated dogs by year of program SD.vaccines.given <- rep(0, times = length(nu.vec)*trials*10*5) dim(SD.vaccines.given) <- c(trials, 5, length(nu.vec), 10) #array of vaccinated dogs by year of program #arrays for transmission throughout time #rows = trials, columns = vacc coverage, dim3 = time SD.store.Sd <- rep(0, times = (length(record.days)*length(nu.vec)*trials)*5) dim(SD.store.Sd) <- c(trials, 5, length(nu.vec), length(record.days)) SD.store.Ed <- rep(0, times = (length(record.days)*length(nu.vec)*trials)*5) dim(SD.store.Ed) <- c(trials, 5, length(nu.vec), length(record.days)) SD.store.Id <- rep(0, times = (length(record.days)*length(nu.vec)*trials)*5) dim(SD.store.Id) <- c(trials, 5, length(nu.vec), length(record.days)) SD.store.Vd <- rep(0, times = (length(record.days)*length(nu.vec)*trials)*5) dim(SD.store.Vd) <- c(trials, 5, length(nu.vec), length(record.days)) SD.store.Sw <- rep(0, times = (length(record.days)*length(nu.vec)*trials)*5) dim(SD.store.Sw) <- c(trials, 5, length(nu.vec), length(record.days)) SD.store.Ew <- rep(0, times = (length(record.days)*length(nu.vec)*trials)*5) dim(SD.store.Ew) <- c(trials, 5, length(nu.vec), length(record.days)) SD.store.Iw <- rep(0, times = (length(record.days)*length(nu.vec)*trials)*5) dim(SD.store.Iw) <- c(trials, 5, length(nu.vec), length(record.days)) SD.store.Nd <- rep(0, times = (length(record.days)*length(nu.vec)*trials)*5) dim(SD.store.Nd) <- c(trials, 5, length(nu.vec), length(record.days)) SD.store.Nw <- rep(0, times = (length(record.days)*length(nu.vec)*trials)*5) dim(SD.store.Nw) <- c(trials, 5, length(nu.vec), length(record.days)) #*************************************** SIMULATIONS ********************************************# for (k in 1:trials) { a <- a.vec[k] b <- b.vec[k] c <- c.vec[k] d <- d.vec[k] sigma <- sigma.vec[k] alpha <- alpha.vec[k] k11 <- SD.k11.vec[k]*alpha/9.38 k12 <- SD.k12.vec[k]*alpha/9.38 k21 <- SD.k21.vec[k]*alpha/Kw k22 <- SD.k22.vec[k]*alpha/Kw gamma <- (a - b)/Kd gammaw <- (c - d)/Kw for (l in 1:5) { for (j in 1:1) { Sd <- rep (0, times = days + start.up) Ed <- rep (0, times = days + start.up) Id <- rep (0, times = days + start.up) Vd <- rep (0, times = days + start.up) Nd <- rep (0, times = days + start.up) Sw <- rep (0, times = days + start.up) Ew <- rep (0, times = days + start.up) Iw <- rep (0, times = days + start.up) Nw <- rep (0, times = days + start.up) rabid.dogs <- rep (0, times = days + start.up) vaccinated.dogs <- rep (0, times = days + start.up) t <- seq (days + start.up) Sd[1] <- 9.5 #susceptible dog density at start Sw[1] <- 2.9 #susceptible wildlife density at start Ed[1] <- 0.025 #exposed dog density at start Ew[1] <- 0.0025 #exposed wildlife density at start Id[1] <- 0.006 #infected dog density at start Iw[1] <- 0.0006 #infected wildlife density at start Vd[1] <- 0 #vaccinated dog density at start Nd[1] <- Sd[1]+Ed[1]+Id[1]+Vd[1] Nw[1] <- Sw[1]+Ew[1]+Iw[1] vaccinated.dogs[1] <- Vd[1] for (i in 2:(days + start.up)) { Sd[i] <- Sd[i-1] + a*(Sd[i-1] + Vd[i-1]) - b*Sd[i-1] - gamma*Sd[i-1]*Nd[i-1] - k11*Id[i-1]*Sd[i-1] - k12*Iw[i-1]*Sd[i-1] Ed[i] <- Ed[i-1] + k11*Id[i-1]*Sd[i-1] + k12*Iw[i-1]*Sd[i-1] - b*Ed[i-1] - sigma*Ed[i-1] - gamma*Ed[i-1]*Nd[i-1] Id[i] <- Id[i-1] + sigma*Ed[i-1] - b*Id[i-1] - alpha*Id[i-1] Vd[i] <- Vd[i-1] - b*Vd[i-1] - gamma*Vd[i-1]*Nd[i-1] Sw[i] <- Sw[i-1] + (c-d)*Sw[i-1] - gammaw*Sw[i-1]*Nw[i-1] - k22*Iw[i-1]*Sw[i-1] - k21*Id[i-1]*Sw[i-1] Ew[i] <- Ew[i-1] + k22*Iw[i-1]*Sw[i-1] + k21*Id[i-1]*Sw[i-1] - d*Ew[i-1] - sigma*Ew[i-1] - gammaw*Ew[i-1]*Nw[i-1] Iw[i] <- Iw[i-1] + sigma*Ew[i-1] - d*Iw[i-1] - alpha*Iw[i-1] Nd[i] <- Sd[i]+Ed[i]+Id[i]+Vd[i] Nw[i] <- Sw[i]+Ew[i]+Iw[i] if (i > start.up + 1) { rabid.dogs[i] <- rabid.dogs[i-1] + sigma*Ed[i-1] } SD.store.Sd[k,l,j,] <- Sd [(start.up + 1):(days + start.up)][record.days] SD.store.Ed[k,l,j,] <- Ed [(start.up + 1):(days + start.up)][record.days] SD.store.Id[k,l,j,] <- Id [(start.up + 1):(days + start.up)][record.days] SD.store.Vd[k,l,j,] <- Vd [(start.up + 1):(days + start.up)][record.days] SD.store.Sw[k,l,j,] <- Sw [(start.up + 1):(days + start.up)][record.days] SD.store.Ew[k,l,j,] <- Ew [(start.up + 1):(days + start.up)][record.days] SD.store.Iw[k,l,j,] <- Iw [(start.up + 1):(days + start.up)][record.days] SD.store.Nd[k,l,j,] <- Nd [(start.up + 1):(days + start.up)][record.days] SD.store.Nw[k,l,j,] <- Nw [(start.up + 1):(days + start.up)][record.days] } SD.rabid [k,l,j,] <- rabid.dogs[year.marks + start.up] Sd.set <- Sd[start.up] Ed.set <- Ed[start.up] Id.set <- Id[start.up] Sw.set <- Sw[start.up] Ew.set <- Ew[start.up] Iw.set <- Iw[start.up] } for(j in 2:length(nu.vec)) { Sd <- rep (0, times = days) Ed <- rep (0, times = days) Id <- rep (0, times = days) Vd <- rep (0, times = days) Nd <- rep (0, times = days) Sw <- rep (0, times = days) Ew <- rep (0, times = days) Iw <- rep (0, times = days) Nw <- rep (0, times = days) rabid.dogs <- rep(0, times = days) vaccinated.dogs <- rep(0, times = days) vaccines.given <- rep(0, times = days) t <- seq (days) nu <- rep (0, times = days) #vaccination vector if (l == 1) nu [round(seq(365/2, 365*9+1, by = 365/2),0)] <- nu.vec[j] if (l == 2) nu [seq(366, 365*9+1, by = 365)] <- nu.vec[j] if (l == 3) nu [seq(365*2+1, 365*9+1, by = 365*2)] <- nu.vec[j] if (l == 4) nu [seq(365*3+1, 365*9+1, by = 365*3)] <- nu.vec[j] Sd[1] <- Sd.set - Sd.set*nu.vec[j] #susceptible dog density at start Sw[1] <- Sw.set #susceptible wildlife density at start Ed[1] <- Ed.set #exposed dog density at start Ew[1] <- Ew.set #exposed wildlife density at start Id[1] <- Id.set #infected dog density at start Iw[1] <- Iw.set #infected wildlife density at start Vd[1] <- Sd.set*nu.vec[j] #vaccinated dog density at start Nd[1] <- Sd[1]+Ed[1]+Id[1]+Vd[1] Nw[1] <- Sw[1]+Ew[1]+Iw[1] vaccinated.dogs[1] <- Vd[1] vaccines.given[1] <- Vd[1] for (i in 2:days) { if (nu[i] > 0) { Sd[i] <- Sd[i-1] + a*(Sd[i-1] + Vd[i-1]) - b*Sd[i-1] - gamma*Sd[i-1]*Nd[i-1] - k11*Id[i-1]*Sd[i-1] - k12*Iw[i-1]*Sd[i-1] - (nu[i]*Nd[i-1] - Vd[i-1]) Ed[i] <- Ed[i-1] + k11*Id[i-1]*Sd[i-1] + k12*Iw[i-1]*Sd[i-1] - b*Ed[i-1] - sigma*Ed[i-1] - gamma*Ed[i-1]*Nd[i-1] Id[i] <- add[i-1]+ Id[i-1] + sigma*Ed[i-1] - b*Id[i-1] - alpha*Id[i-1] Vd[i] <- Vd[i-1] + (nu[i]*Nd[i-1] - Vd[i-1]) - b*Vd[i-1] - gamma*Vd[i-1]*Nd[i-1] vaccinated.dogs[i] <- vaccinated.dogs[i-1] + (nu[i]*Nd[i-1] - Vd[i-1])/((1.03)^(i/365)) #discounted vaccines.given[i] <- vaccines.given[i-1] + (nu[i]*Nd[i-1] - Vd[i-1])/((1.03)^(i/365)) + rv*Vd[i-1]/((1.03)^(i/365)) #discounted }else{ Sd[i] <- Sd[i-1] + a*(Sd[i-1] + Vd[i-1]) - b*Sd[i-1] - gamma*Sd[i-1]*Nd[i-1] - k11*Id[i-1]*Sd[i-1] - k12*Iw[i-1]*Sd[i-1] Ed[i] <- Ed[i-1] + k11*Id[i-1]*Sd[i-1] + k12*Iw[i-1]*Sd[i-1] - b*Ed[i-1] - sigma*Ed[i-1] - gamma*Ed[i-1]*Nd[i-1] Id[i] <- add[i-1]+ Id[i-1] + sigma*Ed[i-1] - b*Id[i-1] - alpha*Id[i-1] Vd[i] <- Vd[i-1] - b*Vd[i-1] - gamma*Vd[i-1]*Nd[i-1] vaccinated.dogs[i] <- vaccinated.dogs[i-1] vaccines.given[i] <- vaccines.given[i-1] } Sw[i] <- Sw[i-1] + (c-d)*Sw[i-1] - gammaw*Sw[i-1]*Nw[i-1] - k22*Iw[i-1]*Sw[i-1] - k21*Id[i-1]*Sw[i-1] Ew[i] <- Ew[i-1] + k22*Iw[i-1]*Sw[i-1] + k21*Id[i-1]*Sw[i-1] - d*Ew[i-1] - sigma*Ew[i-1] - gammaw*Ew[i-1]*Nw[i-1] Iw[i] <- Iw[i-1] + sigma*Ew[i-1] - d*Iw[i-1] - alpha*Iw[i-1] Nd[i] <- Sd[i]+Ed[i]+Id[i]+Vd[i] Nw[i] <- Sw[i]+Ew[i]+Iw[i] rabid.dogs[i] <- rabid.dogs[i-1] + sigma*Ed[i-1]/((1.03)^(i/365)) #discounted SD.store.Sd[k,l,j,] <- Sd[record.days] SD.store.Ed[k,l,j,] <- Ed[record.days] SD.store.Id[k,l,j,] <- Id[record.days] SD.store.Vd[k,l,j,] <- Vd[record.days] SD.store.Sw[k,l,j,] <- Sw[record.days] SD.store.Ew[k,l,j,] <- Ew[record.days] SD.store.Iw[k,l,j,] <- Iw[record.days] SD.store.Nd[k,l,j,] <- Nd[record.days] SD.store.Nw[k,l,j,] <- Nw[record.days] } SD.rabid [k,l,j,] <- rabid.dogs[year.marks] SD.vaccinated [k,l,j,] <- vaccinated.dogs[year.marks] SD.vaccines.given[k,l,j,] <- vaccines.given[year.marks] } } } close(pb) Sys.time() - clock.start #*************************************** SAVE OUTPUT ********************************************# save(SD.rabid, SD.vaccinated, SD.vaccines.given, SD.k11.vec, SD.k12.vec, SD.k21.vec, SD.k22.vec, file = paste0("Serengeti n = 10000", as.character(Sys.Date()),".RData")) #*************************************** COSTS ********************************************# trials <- 10000 vacc.cost.vec <- c(0,41.11,22.08,15.73,12.57,10.69,9.42,8.48,9.79,10.64,14.13) SD.vacc.cost <- rep(0, length(vacc.cost.vec)*trials*5) dim(SD.vacc.cost) <- c(trials, 5, length(vacc.cost.vec)) SD.DALY.gained <- rep(0, length(vacc.cost.vec)*trials*5) dim(SD.DALY.gained) <- c(trials, 5, length(vacc.cost.vec)) dollars.more.spent <- rep(0, length(vacc.cost.vec)*trials*5) dim(dollars.more.spent) <- c(trials, 5, length(vacc.cost.vec)) for (i in 1:trials) { for (j in 1:5) { for (k in 1:11) SD.vacc.cost[i,j,k] <- SD.vaccines.given[i,j,k,10]*vacc.cost.vec[k] } } SD.PEP.cost <- SD.rabid[,,,10]*36.89*2.8635 SD.monetary.cost <- SD.vacc.cost + SD.PEP.cost SD.DALY.cost <- SD.rabid[,,,10]*1.44 SD.rabid.people <- SD.rabid[,,,10]*(.51*(1-.65)*.19) save(SD.rabid.people, file = "SD rabid people.RData") for (i in 1:trials) { for (j in 1:5) { for (k in 1:11) { SD.DALY.gained[i,j,k] <- SD.DALY.cost[i,j,1] - SD.DALY.cost[i,j,k] dollars.more.spent[i,j,k] <- SD.monetary.cost[i,j,k] - SD.monetary.cost[i,j,1] } } } dollars.per.DALY <- dollars.more.spent/SD.DALY.gained #*************************************** EFFICIENT FRONTIER ********************************************# mval <- matrix(10000000000000000000000000000000000, trials, 10) xval <- matrix(0, trials, 10) yval <- matrix(0, trials,10) for (l in 1:trials) { xval[l,] <- arrayInd(which.min(SD.monetary.cost[l,,1:10]), dim(SD.monetary.cost[l,,1:10]))[1] yval[l,] <- arrayInd(which.min(SD.monetary.cost[l,,1:10]), dim(SD.monetary.cost[l,,1:10]))[2] for (k in 1:9) { for (i in 1:4) { for (j in 1:10) { if (SD.monetary.cost[l,i,j] > SD.monetary.cost[l,xval[l,k],yval[l,k]]) { if (SD.DALY.gained[l,xval[l,k],yval[l,k]]-SD.DALY.gained[l,i,j]!=0) { mval2 = (SD.monetary.cost[l,xval[l,k],yval[l,k]]-SD.monetary.cost[l,i,j])/(SD.DALY.gained[l,xval[l,k],yval[l,k]]-SD.DALY.gained[l,i,j]) if (mval2 < mval[l,k+1] & mval2 >= 0) { mval[l,k+1] = mval2 xval[l,k+1] = i yval[l,k+1] = j } } } } } } } for (i in 1:trials) { for (j in 1:9) { x <- xval[i,j] y <- yval[i,j] if (x == 1) { if (y==10) { for (k in (j+1):10) { xval[i,k] <- 0 yval[i,k] <- 0 } } } } } icer <- matrix(0, trials, 10) for (l in 1:trials) { for (i in 10:2) { if (xval[l,i]==0) { icer[l,i] = 0 } else { icer[l,i] = (SD.monetary.cost[l, xval[l, i], yval[l, i]]-SD.monetary.cost[l, xval[l, i-1], yval[l, i-1]])/(SD.DALY.gained[l, xval[l, i], yval[l, i]]-SD.DALY.gained[l, xval[l, i-1], yval[l, i-1]]) } } } #*************************************** NET BENEFITS ********************************************# ceiling.ratio <- c(seq(1, 15000, by = 10)) SD.max.nb <- rep(0, trials*length(ceiling.ratio)) dim(SD.max.nb) <- c(trials, length(ceiling.ratio)) maxX <- matrix(0,trials, length(ceiling.ratio)) maxY <- matrix(0,trials, length(ceiling.ratio)) for (i in 1:length(ceiling.ratio)) { nb <- SD.DALY.gained[,1:4,1:10]-SD.monetary.cost[,1:4,1:10]/ceiling.ratio[i] SD.max.nb[,i] <- apply(nb, 1, which.max) } #*************************************** PROCESS FOR GRAPHS ********************************************# test <- rep(0, times = length(ceiling.ratio)*4*10) dim(test) <- c(length(ceiling.ratio),4,10) for(i in 1:length(ceiling.ratio)) { for(j in 1:trials) { midval <- SD.max.nb[j,i]%%4 if (midval == 0) midval <- 4 test[i, midval, ceiling(SD.max.nb[j,i]/4)] <- test[i, midval, ceiling(SD.max.nb[j,i]/4)] + 1 } } prob <- test/trials SDprob <- prob library(reshape) SDprob <- melt(SDprob) x <- (5/3373)/10*1.5*14036 SDprob$rval <- rep(seq(1, 15000, by = 10), times = 40) freq <- c("6 months", "1 year", "2 years", "3 years") SDprob$freq <- factor(freq[SDprob$X2], levels = freq) percs <- c("0%", "10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%") #percs <- c("10%", "30%", "50%", "70%", "90%") SDprob$txt <- factor(percs[SDprob$X3], levels=percs) SDprob$ind <- paste0(SDprob$freq, SDprob$txt) SDprob$col <- factor(col_pal2[SDprob$freq]) SDprob$shp <- factor(fig[SDprob$txt], levels = fig) save(SDprob, file = paste0("SDprob ", Sys.Date(), ".RData"))