knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
This is an example of how to run a numerical experimental design with landsepi.
library(landsepi)
First create a dataframe containing, for each line, the value of target parameters (i.e. those destined to vary from one simulation to another). The dataframe must also contain columns to store simulation outputs.
For the sake of simplicity, the following numerical design contains only 5 simulations for which 7 target parameters have been arbitrarily chosen to vary and 3 output variable are to be computed. Off course, real numerical designs will include more simulations, in particular those representing factorial experimental plans.
myDesign <- data.frame(nTSpY = c(120, 110, 100, 90, 80) , pI0 = c(0, 1E-4, 5E-4, 1E-3, 1E-2) , infection_rate = c(0.5, 0.4, 0.3, 0.2, 0.1) , id_landscape = 1:5 , aggreg = c(0.07, 0.07, 0.25, 10, 10) , R_efficiency = c(1.00, 0.90, 0.80, 0.70, 0.60) , growth_rate = c(0.1, 0.2, 0.3, 0.1, 0.1) ## create columns to store outputs , durab_MG1 = NA , durab_MG2 = NA , mean_audpc = NA)
Then add an identifier for each simulation and specify a random seed (for stochasticity).
n <- nrow(myDesign) myDesign <- cbind(simul = 1:n, seed = sample(n*100, n), myDesign) myDesign
Create the object simul_params
and a repository to store inputs and outputs of the simulations.
simul_params <- createSimulParams(outputDir = getwd())
Then run the experimental design by updating target parameters with the values stored in myDesign
.
Note that some parameters vary jointly (e.g. landscape and dispersal matrix).
See tutorial on how to run a simple simulation for details on the different functions.
for (i in 1:n){ print(paste("Running simulation", i, "/", n)) ## Set Nyears and nTSpY simul_params <- setTime(simul_params , Nyears = 6 , nTSpY = myDesign$nTSpY[i]) ## update nTSpY ## set seed (to run stochastic replicates) simul_params <- setSeed(simul_params, myDesign$seed[i]) ## update seed ## Pathogen parameters simul_params@ReproSexProb <- logical(0) ## initialize vector of sexual probability (one for each time step) basic_patho_param <- loadPathogen("rust") basic_patho_param$infection_rate <- myDesign$infection_rate[i] ## update inf. rate simul_params <- setPathogen(simul_params, basic_patho_param) ## Initial conditions simul_params <- setInoculum(simul_params, myDesign$pI0[i]) ## update pI0 ## Landscape parameters landscape <- loadLandscape(myDesign$id_landscape[i]) ## update landscape simul_params <- setLandscape(simul_params, landscape) ## Dispersal parameters disp_patho_clonal <- loadDispersalPathogen(myDesign$id_landscape[i])[[1]] ## update dispersal simul_params <- setDispersalPathogen(simul_params, disp_patho_clonal) ## Genes gene1 <- loadGene(name = "MG 1", type = "majorGene") gene1$mutation_prob<-1e-4 gene2 <- loadGene(name = "MG 2", type = "majorGene") gene2$mutation_prob<-1e-4 genes <- data.frame(rbind(gene1, gene2), stringsAsFactors = FALSE) genes$efficiency <- myDesign$R_efficiency[i] ## update resistance efficiency simul_params <- setGenes(simul_params, genes) ## Cultivars cultivar1 <- loadCultivar(name = "Susceptible", type = "growingHost") cultivar2 <- loadCultivar(name = "Resistant1", type = "growingHost") cultivar3 <- loadCultivar(name = "Resistant2", type = "growingHost") cultivars <- data.frame(rbind(cultivar1, cultivar2, cultivar3) , stringsAsFactors = FALSE) cultivars$growth_rate <- myDesign$growth_rate[i] ## update growth rate simul_params <- setCultivars(simul_params, cultivars) ## Allocate genes to cultivars simul_params <- allocateCultivarGenes(simul_params, "Resistant1", c("MG 1")) simul_params <- allocateCultivarGenes(simul_params, "Resistant2", c("MG 2")) ## Allocate cultivars to croptypes croptypes <- loadCroptypes(simul_params, names = c("Susceptible crop" , "Resistant crop 1" , "Resistant crop 2")) croptypes <- allocateCroptypeCultivars(croptypes, "Susceptible crop", "Susceptible") croptypes <- allocateCroptypeCultivars(croptypes, "Resistant crop 1", "Resistant1") croptypes <- allocateCroptypeCultivars(croptypes, "Resistant crop 2", "Resistant2") simul_params <- setCroptypes(simul_params, croptypes) ## Allocate croptypes to landscape rotation_sequence <- croptypes$croptypeID ## No rotation: 1 rotation_sequence element rotation_period <- 0 ## same croptypes every years prop <- c(1/3, 1/3, 1/3) ## croptypes proportions aggreg <- myDesign$aggreg[i] simul_params <- allocateLandscapeCroptypes(simul_params, rotation_period = rotation_period, rotation_sequence = rotation_sequence, rotation_realloc = FALSE, prop = prop, aggreg = aggreg, graphic = FALSE) ## configure outputs outputlist <- loadOutputs(epid_outputs = "audpc", evol_outputs = "durability") simul_params <- setOutputs(simul_params, outputlist) ## Check, (save) and run simulation checkSimulParams(simul_params) # simul_params <- saveDeploymentStrategy(simul_params, overwrite = TRUE) res <- runSimul(simul_params, writeTXT=FALSE, graphic = FALSE) ## Extract outputs myDesign$durab_MG1[i] <- res$evol_outputs$durability[,"MG 1"] myDesign$durab_MG2[i] <- res$evol_outputs$durability[,"MG 2"] myDesign$mean_audpc[i] <- mean(res$epid_outputs$audpc$total) ## Create myDesign.txt at first simulation and then append if (i ==1){ write.table(myDesign[i,], paste(simul_params@OutputDir, "/myDesign.txt", sep="") , append=FALSE, row.names = FALSE, col.names = TRUE) } else { write.table(myDesign[i,], paste(simul_params@OutputDir, "/myDesign.txt", sep="") , append=TRUE, row.names = FALSE, col.names = FALSE) } } myDesign
The argument "keepRawResults" from function runSimul()
allows to keep a set of binary files after the
end of the simulation.
This set of binary files is composed of one file per simulated year and per compartment:
- H: healthy hosts,
- Hjuv: juvenile healthy hosts (for host reproduction),
- L: latently infected hosts,
- I: infectious hosts,
- R: removed hosts,
- P: propagules.
Each file indicates for every time-step the number of individuals in each field, and when appropriate for each host and pathogen genotypes.
## Disable computation of outputs: outputlist <- loadOutputs(epid_outputs = "", evol_outputs = "") simul_params <- setOutputs(simul_params, outputlist) ## Run simulation and keep raw binary files: runSimul(simul_params, writeTXT=FALSE, graphic = FALSE, keepRawResults = TRUE)
The following code loads the content of the binary files in lists (named "H", "Hjuv", "P", "L", "I", "R"). Each element of these lists contains a matrix or an array of the number of individuals associated with each field, host genotype and pathogen genotype for a given time step (i.e. the number of elements of each list is the total number of time steps). Then, users can compute their own output variables using these lists.
## retrieve parameters from the object simul_params path <- simul_params@OutputDir Nyears <- simul_params@TimeParam$Nyears nTSpY <- simul_params@TimeParam$nTSpY nTS <- Nyears * nTSpY ## Total number of time-steps Npoly <- nrow(simul_params@Landscape) Nhost <- nrow(simul_params@Cultivars) Npatho <- prod(simul_params@Genes$Nlevels_aggressiveness) ## Initialise lists H <- as.list(1:nTS) Hjuv <- as.list(1:nTS) P <- as.list(1:nTS) L <- as.list(1:nTS) I <- as.list(1:nTS) R <- as.list(1:nTS) index <- 0 ## Read binary files and store values in the lists as matrices or arrays for (year in 1:Nyears) { binfileH <- file(paste(path, sprintf("/H-%02d", year), ".bin", sep = ""), "rb") H.tmp <- readBin(con = binfileH, what = "int", n = Npoly * Nhost * nTSpY, size = 4 , signed = T, endian = "little") close(binfileH) binfileHjuv = file(paste(path, sprintf("/Hjuv-%02d", year), ".bin",sep=""), "rb") Hjuv.tmp <- readBin(con=binfileHjuv, what="int", n=Npoly*Nhost*nTSpY, size = 4 , signed=T,endian="little") close(binfileHjuv) binfileP <- file(paste(path, sprintf("/P-%02d", year), ".bin", sep = ""), "rb") P.tmp <- readBin(con = binfileP, what = "int", n = Npoly * Npatho * nTSpY, size = 4 , signed = T, endian = "little") close(binfileP) binfileL <- file(paste(path, sprintf("/L-%02d", year), ".bin", sep = ""), "rb") L.tmp <- readBin(con = binfileL, what = "int", n = Npoly * Npatho * Nhost * nTSpY , size = 4 , signed = T, endian = "little") close(binfileL) binfileI <- file(paste(path, sprintf("/I-%02d", year), ".bin", sep = ""), "rb") I.tmp <- readBin(con = binfileI, what = "int", n = Npoly * Npatho * Nhost * nTSpY , size = 4 , signed = T, endian = "little") close(binfileI) binfileR <- file(paste(path, sprintf("/R-%02d", year), ".bin", sep = ""), "rb") R.tmp <- readBin(con = binfileR, what = "int", n = Npoly * Npatho * Nhost * nTSpY , size = 4 , signed = T, endian = "little") close(binfileR) ## Convert vectors in matrices or arrays for (t in 1:nTSpY) { H[[t + index]] <- matrix(H.tmp[((Nhost * Npoly) * (t-1)+1):(t * (Nhost * Npoly))] , ncol = Nhost, byrow = T) Hjuv[[t + index]] <- matrix(Hjuv.tmp[((Nhost*Npoly)*(t-1)+1):(t*(Nhost*Npoly))] , ncol=Nhost,byrow=T) P[[t + index]] <- matrix(P.tmp[((Npatho * Npoly) * (t-1)+1):(t * (Npatho * Npoly))] , ncol = Npatho, byrow = T) L[[t + index]] <- array(data = L.tmp[((Npatho * Npoly * Nhost) * (t-1)+1):(t * (Npatho * Npoly * Nhost))] , dim = c(Nhost, Npatho, Npoly)) I[[t + index]] <- array(data = I.tmp[((Npatho * Npoly * Nhost) * (t-1)+1):(t * (Npatho * Npoly * Nhost))] , dim = c(Nhost, Npatho, Npoly)) R[[t + index]] <- array(data = R.tmp[((Npatho * Npoly * Nhost) * (t-1)+1):(t * (Npatho * Npoly * Nhost))] , dim = c(Nhost, Npatho, Npoly)) } index <- index + nTSpY }
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.