2 - Running a numerical experimental design

This is an example of how to run a numerical experimental design with landsepi.

library(landsepi)

Write the experimental design

First create a dataframe containing, for each line, the value of target parameters (i.e. those destined to vary from one simulation to anoter). 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
#>   simul seed nTSpY   pI0 infection_rate id_landscape aggreg R_efficiency
#> 1     1  472   120 0e+00            0.5            1   0.07          1.0
#> 2     2  334   110 1e-04            0.4            2   0.07          0.9
#> 3     3   71   100 5e-04            0.3            3   0.25          0.8
#> 4     4  384    90 1e-03            0.2            4  10.00          0.7
#> 5     5  315    80 1e-02            0.1            5  10.00          0.6
#>   growth_rate durab_MG1 durab_MG2 mean_audpc
#> 1         0.1        NA        NA         NA
#> 2         0.2        NA        NA         NA
#> 3         0.3        NA        NA         NA
#> 4         0.1        NA        NA         NA
#> 5         0.1        NA        NA         NA

Run the simulations

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
    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 <- loadDispersalPathogen(myDesign$id_landscape[i])  ## update dispersal
    simul_params <- setDispersalPathogen(simul_params, disp_patho)
    
    ## Genes
    gene1 <- loadGene(name = "MG 1", type = "majorGene")
    gene2 <- loadGene(name = "MG 2", type = "majorGene")
    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

Compute your own output variables

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,
- 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 <- 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
}