Skip to contents

This document is still work in progress

OpenMalaria is a powerful tool for simulating and studying the epidemiology of malaria. This R package allows the use of OpenMalaria from R.

Prerequisites

Make sure that SQLite and OpenMalaria are both installed and available on the system. Please note: Even though the package should work on Windows, it is developed and tested under Linux.

To install the package, run

devtools::install_github("SwissTPH/r-openMalariaUtilities")

Introduction

Within this tutorial, we will create and run example scenarios. Each scenario needs to be defined in a XML file which serves as an input for OpenMalaria (OM). OpenMalariaUtilities (OMU) contains many functions which help to generate these input files. In this tutorial we will re-create the example from OpenMalaria itself.

XML file generation

A large part of OMU is dedicated at the generation of valid OM XML files. Before we start with the generation of these files, we will explain the folder layout of your projects.

In the folder hierarchy, your working or project directory is the top level. Within this folder, you can have one or more folders which correspond to the names of experiments. The simulations of these experiments are controlled via one or multiple R scripts, placed in the project directory. Please note: Most of the directories can be modified via function arguments, but for the sake of simplicity we using the defaults here.

MyProject
├── experiment_1
├── experiment_2
├── ...
├── script_1.R
├── script_2.R
└── ...

Each of the experiment folders will have a similar layout and content.

MyProject/experiment_1
├── cache/
├── logs/
├── outputs/
├── scenarios/
├── autoRegressionParameters.csv
├── densities.csv
├── base.xml
└── scenario_44.xsd

You do not need to worry about the creation of the folders and files. OMU will create necessary folders and download required files for you.

The cache folder contains R objects which are used by OMU to store reusable information. This will speed up certain operations. Log files are stored in the logs folder (which will contain sub-folders for each step in the process) and OpenMalaria’s raw output will be outputs. The XML files for each simulation scenario are stored in the scenarios folder.

Starting the project

In your R script, you will probably start with library(openMalariaUtilities). Thus, the beginning could look like

library(openMalariaUtilities)

## Basic skeleton
baseList <- list(
  ## Mandatory
  expName = "example",
  ## Mandatory
  OMVersion = 44L,
  ## Optional
  ## analysisNo = 1L,
  ## Mandatory
  demography = list(),
  monitoring = list(),
  interventions = list(),
  healthSystem = list(),
  entomology = list(),
  ## These are optional for OM
  ## parasiteGenetics = list(),
  ## pharmacology = list(),
  ## diagnostics = list(),
  model = list()
)

For XML file generation, OMU uses a list as input (here: baseList). This list has a very similar structure as the XML file required by OM but is easier to edit from R compared to XML. The above snippet corresponds to the minimal skeleton of this input list. expName is the name of your experiment and OMVersion the version of OpenMalaria. Currently, only version 44 is supported. Also note, that the OM version needs to be an integer.

Each of the six nested lists demography, monitoring, interventions, healthSystem, entomology and model need to be filled for OM.

Demography

There are two ways to add the data to the demography entry:

  1. Add it manually, e.g. baseList[["demography"]] <- list(name = "Ifakara", ...)
  2. Use the provided defineDemography() function

We will use the second option in this tutorial. You can find various defineXY functions which are each designed to add a specific entry to the list. However, you always have the option to write parts or change specific entries manually (See also extractList()).

demo_data <- data.frame(
  poppercent = c(
    3.474714994, 12.76004028, 14.52151394, 12.75565434, 10.83632374,
    8.393312454, 7.001421452, 5.800587654, 5.102136612, 4.182561874,
    3.339409351, 2.986112356, 2.555766582, 2.332763433, 1.77400255,
    1.008525491, 0.74167341, 0.271863401, 0.161614642
  ),
  upperbound = c(
    1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90
  )
)


baseList <- defineDemography(
  baseList,
  name = "Ifakara",
  popSize = 4000L,
  maximumAgeYrs = 90,
  lowerbound = 0,
  poppercent = demo_data$poppercent,
  upperbound = demo_data$upperbound
)

str(baseList, max.level = 2)
#> List of 8
#>  $ expName      : chr "example"
#>  $ OMVersion    : int 44
#>  $ demography   :List of 4
#>   ..$ name         : chr "Ifakara"
#>   ..$ popSize      : int 4000
#>   ..$ maximumAgeYrs: num 90
#>   ..$ ageGroup     :List of 20
#>  $ monitoring   : list()
#>  $ interventions: list()
#>  $ healthSystem : list()
#>  $ entomology   : list()
#>  $ model        : list()

In the above snippet, we started by creating a data frame which contains the demographic data. This in turn was used as input for defineDemography() among other parameters. The function will take care that the information is added at the correct position and also check, that the correct data types are used.

Monitoring

The monitoring section controls the output of OpenMalaria, which is generated during the conducted surveys.

baseList[["monitoring"]] <- list(
  name = "Quarterly Surveys",
  ## Mandatory, different from OM schema
  startDate = "2000-01-01",
  continuous = monitoringContinuousGen(
    period = 1,
    options = list(
      name = c("input EIR", "simulated EIR", "human infectiousness", "N_v0",
               "immunity h", "immunity Y", "new infections",
               "num transmitting humans", "ITN coverage", "GVI coverage", "alpha",
               "P_B", "P_C*P_D"),
      value = c("true", "true", "true", "true", "true", "true", "true", "true",
                "true", "true", "true", "false", "false")
    )
  ),
  SurveyOptions = monitoringSurveyOptionsGen(
    options = list(
      name = c("nHost", "nPatent", "nUncomp", "nSevere", "nDirDeaths",
               "inputEIR", "simulatedEIR"),
      value = c("true", "true", "true", "true", "false", "true", "true")
    )
  )
)

We started the monitoring definition by giving it a name and defining a start date. Here, this is mandatory because it allows us to use dates as input for OpenMalaria. It defines the starting point of the survey period.

For the options, two generator functions are used which take the name of the option and its value as input.

baseList[["monitoring"]] <- c(
  baseList[["monitoring"]],
  list(
    surveys =  monitoringSurveyTimesGen(
      detectionLimit = 40, startDate = "2000-01-01", endDate = "2020-01-01",
      interval = "quarterly", simStart = "1980-01-01"
    ),
    ageGroup = surveyAgeGroupsGen(
      lowerbound = 0,
      upperbounds = c(5, 90)
    )
  )
)

str(baseList[["monitoring"]], max.level = 1)
#> List of 6
#>  $ name         : chr "Quarterly Surveys"
#>  $ startDate    : chr "2000-01-01"
#>  $ continuous   :List of 14
#>  $ SurveyOptions:List of 7
#>  $ surveys      :List of 86
#>  $ ageGroup     :List of 3

Similar to the age groups in the demography section, surveyAgeGroupsGen() adds the desired age groups for the surveys.

The definition of the survey times itself is done above. Please note, that you really should use monitoringSurveyTimesGen() as this function will make sure your desired survey points are valid (and adjust them if necessary) and will store them for later use in the cache.

We can inspect the cache to see if our survey times have been generated correctly.

listCache()
#> 
#> Name: mon_ageGroups  Timestamp: 2024-09-06 12:52:36 
#> Class: list 
#> Value:
#>  $lowerbound
#>  [1] 0
#>  
#>  $upperbounds
#>  [1]  5 90
#>  
#>  
#>  -------------------------------------------------------------------------------- 
#> 
#> Name: simStart   Timestamp: 2024-09-06 12:52:36 
#> Class: character 
#> Value:
#>  [1] "1980-01-01"
#>  
#>  -------------------------------------------------------------------------------- 
#> 
#> Name: surveyTimes    Timestamp: 2024-09-06 12:52:36 
#> Class: data.table data.frame 
#> Value:
#>      number       date   day timestep
#>       <int>     <Date> <num>    <num>
#>   1:      1 1999-12-27  7295     1459
#>   2:      2 2000-03-27  7385     1477
#>   3:      3 2000-06-30  7480     1496
#>   4:      4 2000-09-28  7570     1514
#>   5:      5 2000-12-27  7660     1532
#>   6:      6 2001-03-27  7750     1550
#>   7:      7 2001-06-30  7845     1569
#>   8:      8 2001-09-28  7935     1587
#>   9:      9 2001-12-27  8025     1605
#>  10:     10 2002-03-27  8115     1623
#>  11:     11 2002-06-30  8210     1642
#>  12:     12 2002-09-28  8300     1660
#>  13:     13 2002-12-27  8390     1678
#>  14:     14 2003-03-27  8480     1696
#>  15:     15 2003-06-30  8575     1715
#>  16:     16 2003-09-28  8665     1733
#>  17:     17 2003-12-27  8755     1751
#>  18:     18 2004-03-27  8845     1769
#>  19:     19 2004-06-30  8940     1788
#>  20:     20 2004-09-28  9030     1806
#>  21:     21 2004-12-27  9120     1824
#>  22:     22 2005-03-27  9210     1842
#>  23:     23 2005-06-30  9305     1861
#>  24:     24 2005-09-28  9395     1879
#>  25:     25 2005-12-27  9485     1897
#>  26:     26 2006-03-27  9575     1915
#>  27:     27 2006-06-30  9670     1934
#>  28:     28 2006-09-28  9760     1952
#>  29:     29 2006-12-27  9850     1970
#>  30:     30 2007-03-27  9940     1988
#>  31:     31 2007-06-30 10035     2007
#>  32:     32 2007-09-28 10125     2025
#>  33:     33 2007-12-27 10215     2043
#>  34:     34 2008-03-27 10305     2061
#>  35:     35 2008-06-30 10400     2080
#>  36:     36 2008-09-28 10490     2098
#>  37:     37 2008-12-27 10580     2116
#>  38:     38 2009-03-27 10670     2134
#>  39:     39 2009-06-30 10765     2153
#>  40:     40 2009-09-28 10855     2171
#>  41:     41 2009-12-27 10945     2189
#>  42:     42 2010-03-27 11035     2207
#>  43:     43 2010-06-30 11130     2226
#>  44:     44 2010-09-28 11220     2244
#>  45:     45 2010-12-27 11310     2262
#>  46:     46 2011-03-27 11400     2280
#>  47:     47 2011-06-30 11495     2299
#>  48:     48 2011-09-28 11585     2317
#>  49:     49 2011-12-27 11675     2335
#>  50:     50 2012-03-27 11765     2353
#>  51:     51 2012-06-30 11860     2372
#>  52:     52 2012-09-28 11950     2390
#>  53:     53 2012-12-27 12040     2408
#>  54:     54 2013-03-27 12130     2426
#>  55:     55 2013-06-30 12225     2445
#>  56:     56 2013-09-28 12315     2463
#>  57:     57 2013-12-27 12405     2481
#>  58:     58 2014-03-27 12495     2499
#>  59:     59 2014-06-30 12590     2518
#>  60:     60 2014-09-28 12680     2536
#>  61:     61 2014-12-27 12770     2554
#>  62:     62 2015-03-27 12860     2572
#>  63:     63 2015-06-30 12955     2591
#>  64:     64 2015-09-28 13045     2609
#>  65:     65 2015-12-27 13135     2627
#>  66:     66 2016-03-27 13225     2645
#>  67:     67 2016-06-30 13320     2664
#>  68:     68 2016-09-28 13410     2682
#>  69:     69 2016-12-27 13500     2700
#>  70:     70 2017-03-27 13590     2718
#>  71:     71 2017-06-30 13685     2737
#>  72:     72 2017-09-28 13775     2755
#>  73:     73 2017-12-27 13865     2773
#>  74:     74 2018-03-27 13955     2791
#>  75:     75 2018-06-30 14050     2810
#>  76:     76 2018-09-28 14140     2828
#>  77:     77 2018-12-27 14230     2846
#>  78:     78 2019-03-27 14320     2864
#>  79:     79 2019-06-30 14415     2883
#>  80:     80 2019-09-28 14505     2901
#>  81:     81 2019-12-27 14595     2919
#>  82:     82 2020-03-27 14685     2937
#>  83:     83 2020-06-30 14780     2956
#>  84:     84 2020-09-28 14870     2974
#>  85:     85 2020-12-27 14960     2992
#>      number       date   day timestep
#>  
#>  --------------------------------------------------------------------------------

We can notice two things:

  1. We are expecting 81 surveys, however we have 85
  2. The dates have been adjusted

OMU adds one more year (here: 4 surveys) to the survey period. This is done in order to ensure that all deployments show their effect. The dates are adjusted so they match OM’s internal time steps. For more details on this, see here.

Furthermore, the number of days or time steps is counted from the simulation start simStart, which is here the 1st January 1980. It is recommended to start the simulation 20 - 50 years before the first survey.

Interventions

In this section we will define actual measures which will be used throughout the simulation.

baseList[["interventions"]] <- list(
  name = "test",
  human = list(
    component = list(
      id = "GVI",
      name = "DDT test",
      GVI = list(
        decay = list(
          L = "0.5",
          "function" = "exponential"
        ),
        anophelesParams = list(
          mosquito = "gambiae_ss", propActive = "1",
          deterrency = list(value = "0.56"),
          preprandialKillingEffect = list(value = "0"),
          postprandialKillingEffect = list(value = "0.24")
        )
      )
    ),
    deployment = list(
      name = "DDT test",
      component = list(id = "GVI"),
      timed = list(
        deploy = list(coverage = "0.9", time = "10y")
      )
    )
  )
)

Health system

baseList[["healthSystem"]] <- list(
  ImmediateOutcomes = list(
    name = "Tanzania ACT",
    drugRegimen = list(
      firstLine = "ACT",
      inpatient = "QN",
      secondLine = "ACT"
    ),
    initialACR = list(
      ACT = list(value = 0.85),
      QN = list(value = 0.998),
      selfTreatment = list(value = 0.63)
    ),
    compliance = list(
      ACT = list(value = 0.9),
      selfTreatment = list(value = 0.85)
    ),
    nonCompliersEffective = list(
      ACT = list(value = 0),
      selfTreatment = list(value = 0)
    ),
    treatmentActions = list(
      ACT = list(
        name = "clear blood-stage infections",
        clearInfections = list(
          stage = "blood",
          timesteps = "1"
        )
      ),
      QN = list(
        name = "clear blood-stage infections",
        clearInfections = list(
          stage = "blood",
          timesteps = "1"
        )
      )
    ),
    pSeekOfficialCareUncomplicated1 = list(value = 0.04),
    pSelfTreatUncomplicated = list(value = 0.01),
    pSeekOfficialCareUncomplicated2 = list(value = 0.04),
    pSeekOfficialCareSevere = list(value = 0.48)
  ),
  CFR = list(
    group = list(lowerbound = 0, value = 0.09189),
    group = list(lowerbound = 0.25, value = 0.0810811),
    group = list(lowerbound = 0.75, value = 0.0648649),
    group = list(lowerbound = 1.5, value = 0.0689189),
    group = list(lowerbound = 2.5, value = 0.0675676),
    group = list(lowerbound = 3.5, value = 0.0297297),
    group = list(lowerbound = 4.5, value = 0.0459459),
    group = list(lowerbound = 7.5, value = 0.0945946),
    group = list(lowerbound = 12.5, value = 0.1243243),
    group = list(lowerbound = 15, value = 0.1378378)
  ),
  ## Mandatory
  pSequelaeInpatient = list(
    interpolation = "none",
    group = list(lowerbound = "0.0", value = 0.0132),
    group = list(lowerbound = "5.0", value = 0.005)
  )
)

Entomology

baseList[["entomology"]] <- list(
  mode = "dynamic", name = "Namawala",
  vector = list(
    anopheles = list(
      mosquito = "gambiae_ss", propInfected = 0.078, propInfectious = "0.021",
      seasonality = list(
        annualEIR = "16", input = "EIR",
        fourierSeries = list(
          EIRRotateAngle = "0",
          coeffic = list(
            a = "0.8968", b = "2.678"
          ),
          coeffic = list(a = "-0.4551", b = "2.599")
        )
      ),
      mosq = list(
        minInfectedThreshold = "0.001",
        mosqRestDuration = list(value = "3"),
        extrinsicIncubationPeriod = list(value = "11"),
        mosqLaidEggsSameDayProportion = list(value = "0.313"),
        mosqSeekingDuration = list(value = "0.33"),
        mosqSurvivalFeedingCycleProbability = list(value = "0.623"),
        availability = list(distr = "const"),
        mosqProbBiting = list(mean = "0.95", variance = "0"),
        mosqProbFindRestSite = list(mean = "0.95", variance = "0"),
        mosqProbResting = list(mean = "0.99", variance = "0"),
        mosqProbOvipositing = list(value = "0.88"),
        mosqHumanBloodIndex = list(value = "0.939")
      ),
      nonHumanHosts = list(
        name = "unprotectedAnimals",
        mosqRelativeEntoAvailability = list(value = "1.0"),
        mosqProbBiting = list(value = "0.95"),
        mosqProbFindRestSite = list(value = "0.95"),
        mosqProbResting = list(value = "0.99")
      )
    ),
    nonHumanHosts = list(name = "unprotectedAnimals", number = "1.0")
  )
)

Model

baseList[["model"]] <- list(
  ModelOptions = list(
    option = list(name = "LOGNORMAL_MASS_ACTION", value = "true"),
    option = list(name = "NO_PRE_ERYTHROCYTIC", value = "true"),
    option = list(name = "INNATE_MAX_DENS", value = "false"),
    option = list(name = "INDIRECT_MORTALITY_FIX", value = "false"),
    option = list(name = "NON_MALARIA_FEVERS", value = "true")
  ),
  clinical = list(
    healthSystemMemory = "6",
    NonMalariaFevers = list(
      incidence = list(
        group = list(lowerbound = "0", value = "0.322769924518357"),
        group = list(lowerbound = "0", value = "0.308520194304172"),
        group = list(lowerbound = "1", value = "0.279441774808493"),
        group = list(lowerbound = "2", value = "0.250431781111273"),
        group = list(lowerbound = "3", value = "0.223285859756841"),
        group = list(lowerbound = "4", value = "0.199298352451799"),
        group = list(lowerbound = "5", value = "0.179376872365614"),
        group = list(lowerbound = "6", value = "0.163623659390782"),
        group = list(lowerbound = "7", value = "0.152227726923469"),
        group = list(lowerbound = "8", value = "0.145022785567758"),
        group = list(lowerbound = "9", value = "0.141493087461765"),
        group = list(lowerbound = "10", value = "0.140473293219353"),
        group = list(lowerbound = "11", value = "0.141109775159515"),
        group = list(lowerbound = "12", value = "0.142644475217328"),
        group = list(lowerbound = "13", value = "0.144335079395766"),
        group = list(lowerbound = "14", value = "0.145964032924869"),
        group = list(lowerbound = "15", value = "0.147708915135714"),
        group = list(lowerbound = "16", value = "0.149731543445568"),
        group = list(lowerbound = "17", value = "0.151887428568276"),
        group = list(lowerbound = "18", value = "0.154060663485195"),
        group = list(lowerbound = "19", value = "0.156179169710494"),
        group = list(lowerbound = "20", value = "0.158135015380583"),
        group = list(lowerbound = "21", value = "0.159704766482219"),
        group = list(lowerbound = "22", value = "0.160807788387655"),
        group = list(lowerbound = "23", value = "0.161427976448279"),
        group = list(lowerbound = "24", value = "0.161620429119137"),
        group = list(lowerbound = "25", value = "0.16144021875986"),
        group = list(lowerbound = "26", value = "0.160943264630612"),
        group = list(lowerbound = "27", value = "0.160217573697398"),
        group = list(lowerbound = "28", value = "0.159422614374451"),
        group = list(lowerbound = "29", value = "0.158542519631641"),
        group = list(lowerbound = "30", value = "0.157501217628248"),
        group = list(lowerbound = "31", value = "0.156175160594841"),
        group = list(lowerbound = "32", value = "0.154402302191411"),
        group = list(lowerbound = "33", value = "0.152102040636481"),
        group = list(lowerbound = "34", value = "0.14921450014676"),
        group = list(lowerbound = "35", value = "0.145714433541659"),
        group = list(lowerbound = "36", value = "0.141800502067518"),
        group = list(lowerbound = "37", value = "0.137916853907569"),
        group = list(lowerbound = "38", value = "0.134503529382102"),
        group = list(lowerbound = "39", value = "0.131746276580642"),
        group = list(lowerbound = "40", value = "0.12969902537497"),
        group = list(lowerbound = "41", value = "0.128398077347679"),
        group = list(lowerbound = "42", value = "0.127864136551891"),
        group = list(lowerbound = "43", value = "0.12804497197004"),
        group = list(lowerbound = "44", value = "0.128894055047661"),
        group = list(lowerbound = "45", value = "0.130350838992718"),
        group = list(lowerbound = "46", value = "0.132286605622701"),
        group = list(lowerbound = "47", value = "0.134599921072495"),
        group = list(lowerbound = "48", value = "0.137212726976988"),
        group = list(lowerbound = "49", value = "0.140035253913284"),
        group = list(lowerbound = "50", value = "0.142934573453621"),
        group = list(lowerbound = "51", value = "0.145830221511879"),
        group = list(lowerbound = "52", value = "0.148674810561069"),
        group = list(lowerbound = "53", value = "0.151497963594518"),
        group = list(lowerbound = "54", value = "0.15438856687865"),
        group = list(lowerbound = "55", value = "0.157403790093505"),
        group = list(lowerbound = "56", value = "0.16059513222516"),
        group = list(lowerbound = "57", value = "0.16402433342886"),
        group = list(lowerbound = "58", value = "0.16770481415944"),
        group = list(lowerbound = "59", value = "0.171626873047865"),
        group = list(lowerbound = "60", value = "0.175748327054247"),
        group = list(lowerbound = "61", value = "0.180030857856799"),
        group = list(lowerbound = "62", value = "0.184411365583771"),
        group = list(lowerbound = "63", value = "0.188816421789366"),
        group = list(lowerbound = "64", value = "0.19316997803338"),
        group = list(lowerbound = "65", value = "0.197435603275487"),
        group = list(lowerbound = "66", value = "0.201578808813379"),
        group = list(lowerbound = "67", value = "0.205556806881398"),
        group = list(lowerbound = "68", value = "0.209307183457343"),
        group = list(lowerbound = "69", value = "0.212783260344084"),
        group = list(lowerbound = "70", value = "0.215944154621391"),
        group = list(lowerbound = "71", value = "0.218749275266548"),
        group = list(lowerbound = "72", value = "0.221187990639016"),
        group = list(lowerbound = "73", value = "0.223361260399378"),
        group = list(lowerbound = "74", value = "0.225363436789592"),
        group = list(lowerbound = "75", value = "0.227254280093211"),
        group = list(lowerbound = "76", value = "0.229084576349576"),
        group = list(lowerbound = "77", value = "0.230891971097789"),
        group = list(lowerbound = "78", value = "0.232690225166173"),
        group = list(lowerbound = "79", value = "0.234484973338876"),
        group = list(lowerbound = "80", value = "0.236276361586796"),
        group = list(lowerbound = "81", value = "0.238064394629696"),
        group = list(lowerbound = "82", value = "0.239849077182917"),
        group = list(lowerbound = "83", value = "0.241630413957381"),
        group = list(lowerbound = "84", value = "0.243408409659591"),
        group = list(lowerbound = "85", value = "0.245183068991633"),
        group = list(lowerbound = "86", value = "0.246954396651183"),
        group = list(lowerbound = "87", value = "0.248722397331501"),
        group = list(lowerbound = "88", value = "0.250487075721441"),
        group = list(lowerbound = "89", value = "0.252248436505447"),
        group = list(lowerbound = "90", value = "0.253127874257909")
      )
    )
  ),
  human = list(
    availabilityToMosquitoes = list(
      group = list(lowerbound = "0", value = "0.7076"),
      group = list(lowerbound = "0", value = "0.8538"),
      group = list(lowerbound = "5", value = "1.0"),
      group = list(lowerbound = "5", value = "1.0")
    )
  ),
  parameters = list(
    interval = "5", iseed = "1", latentp = "3",
    parameter = list(include = "false", name = "'-ln(1-Sinf)'", number = "1", value = "0.050736"),
    parameter = list(include = "false", name = "Estar", number = "2", value = "0.03247"),
    parameter = list(include = "false", name = "Simm", number = "3", value = "0.138161050830301"),
    parameter = list(include = "false", name = "Xstar_p", number = "4", value = "1514.385853233699891"),
    parameter = list(include = "false", name = "gamma_p", number = "5", value = "2.03692533424484"),
    parameter = list(include = "false", name = "sigma2i", number = "6", value = "10.173598698525799"),
    parameter = list(include = "false", name = "CumulativeYstar", number = "7", value = "35158523.31132510304451"),
    parameter = list(include = "false", name = "CumulativeHstar", number = "8", value = "97.334652723897705"),
    parameter = list(include = "false", name = "'-ln(1-alpha_m)'", number = "9", value = "2.33031045876193"),
    parameter = list(include = "false", name = "decay_m", number = "10", value = "2.53106547375805"),
    parameter = list(include = "false", name = "sigma2_0", number = "11", value = "0.655747311168152"),
    parameter = list(include = "false", name = "Xstar_v", number = "12", value = "0.916181104713054"),
    parameter = list(include = "false", name = "Ystar2", number = "13", value = "6502.26335600001039"),
    parameter = list(include = "false", name = "alpha", number = "14", value = "142601.912520000012591"),
    parameter = list(include = "false", name = "Density bias (non Garki)", number = "15", value = "0.177378570987455"),
    parameter = list(include = "false", name = "        sigma2        ", number = "16", value = "0.05"),
    parameter = list(include = "false", name = "log oddsr CF community", number = "17", value = "0.736202"),
    parameter = list(include = "false", name = "Indirect risk cofactor", number = "18", value = "0.018777338"),
    parameter = list(include = "false", name = "Non-malaria infant mortality", number = "19", value = "49.539046599999999"),
    parameter = list(include = "false", name = "Density bias (Garki)", number = "20", value = "4.79610772546704"),
    parameter = list(include = "false", name = "Severe Malaria Threshhold", number = "21", value = "784455.599999999976717"),
    parameter = list(include = "false", name = "Immunity Penalty", number = "22", value = "1"),
    parameter = list(include = "false", name = "Immune effector decay", number = "23", value = "0"),
    parameter = list(include = "false", name = "comorbidity intercept", number = "24", value = "0.0968"),
    parameter = list(include = "false", name = "Ystar half life", number = "25", value = "0.275437402"),
    parameter = list(include = "false", name = "Ystar1", number = "26", value = "0.596539864"),
    parameter = list(include = "false", name = "Asexual immunity decay", number = "27", value = "0"),
    parameter = list(include = "false", name = "Ystar0", number = "28", value = "296.302437899999973"),
    parameter = list(include = "false", name = "Idete multiplier", number = "29", value = "2.797523626"),
    parameter = list(include = "false", name = "critical age for comorbidity", number = "30", value = "0.117383")
  )
)

Prepare to run the simulations

## Setup dirs
setupDirs(experimentName = "exp_test", replace = TRUE)

## Copy necessary OpenMalaria files.
setupOM()
#> [1] TRUE

setupDirs() creates the experiment folder given as an argument in the project directory. This also includes all the sub-folders which were mentioned in the introduction .

We can verify this by checking the layout of the experiment directory.

list.dirs("exp_test")
#> [1] "exp_test"                 "exp_test/cache"          
#> [3] "exp_test/logs"            "exp_test/logs/results"   
#> [5] "exp_test/logs/scenarios"  "exp_test/logs/simulation"
#> [7] "exp_test/outputs"         "exp_test/scenarios"

setupOM() will make sure that all required files for OpenMalaria to run are in place. This includes the correct .xsd schema file (according to the specified OM version) as well as further supplementary files.

list.files("exp_test")
#> [1] "autoRegressionParameters.csv" "cache"                       
#> [3] "densities.csv"                "logs"                        
#> [5] "outputs"                      "scenario_44.xsd"             
#> [7] "scenarios"

Create XML file

As the name implies, createBaseXml() will translate our input list to a XML document and write it to disk.

## Create base XML file
createBaseXml(baseList, replace = TRUE)
#> Warning in createBaseXml(baseList, replace = TRUE): startDate 2000-01-01 is not
#> the same as simStart 1980-01-01. This can cause survey times to not appear!
readLines(file("exp_test/exp_test_base.xml", "r"), n = 10)
#>  [1] "<?xml version=\"1.0\" encoding=\"UTF-8\"?>"                                                                                                                                                                                                          
#>  [2] "<om:scenario xmlns:om=\"http://openmalaria.org/schema/scenario_44\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" name=\"exp_test\" schemaVersion=\"44\" xsi:schemaLocation=\"http://openmalaria.org/schema/scenario_44 scenario_44.xsd\">"
#>  [3] "  <demography name=\"Ifakara\" popSize=\"4000\" maximumAgeYrs=\"90\">"                                                                                                                                                                               
#>  [4] "    <ageGroup lowerbound=\"0\">"                                                                                                                                                                                                                     
#>  [5] "      <group poppercent=\"3.474714994\" upperbound=\"1\"/>"                                                                                                                                                                                          
#>  [6] "      <group poppercent=\"12.76004028\" upperbound=\"5\"/>"                                                                                                                                                                                          
#>  [7] "      <group poppercent=\"14.52151394\" upperbound=\"10\"/>"                                                                                                                                                                                         
#>  [8] "      <group poppercent=\"12.75565434\" upperbound=\"15\"/>"                                                                                                                                                                                         
#>  [9] "      <group poppercent=\"10.83632374\" upperbound=\"20\"/>"                                                                                                                                                                                         
#> [10] "      <group poppercent=\"8.393312454\" upperbound=\"25\"/>"

Run simulations

In order to run the simulation, we can simply call runSimulations().

## Run simulations
runSimulations()
#> [1] "Running scenario [1/1]"

By default, the output in the R console is fairly minimal. More detailed logs can be found in the logs/simulation folder.

list.files("exp_test/logs/simulation")
#> [1] "exp_test_base.log"

OpenMalaria generates two output files which are both placed in the outputs folder:

  1. The *_cts.txt file
read.table("exp_test/outputs/exp_test_base_cts.txt", header = TRUE, sep = "\t", skip = 1, nrows = 10)
#>    timestep  input.EIR simulated.EIR human.infectiousness N_v0.gambiae_ss.
#> 1         0 0.00995606    0.00995049            0.0179202          30147.4
#> 2         1 0.01958060    0.15390600            0.0175667          39181.0
#> 3         2 0.03850920    0.24658900            0.0171114          46368.5
#> 4         3 0.07460740    0.37191900            0.0167741          49838.5
#> 5         4 0.14030700    0.52390000            0.0164523          48641.0
#> 6         5 0.25256000    0.68747000            0.0165751          43197.7
#> 7         6 0.42953400    0.84148200            0.0188351          35064.0
#> 8         7 0.68226100    0.96738200            0.0231314          26185.7
#> 9         8 1.00214000    1.07333000            0.0296558          18146.4
#> 10        9 1.35047000    1.18202000            0.0377866          11790.1
#>    immunity.h immunity.Y new.infections num.transmitting.humans ITN.coverage
#> 1     256.846   11899700              0                    2613            0
#> 2     256.817   11900500            400                    2590            0
#> 3     256.856   11895200            692                    2572            0
#> 4     256.934   11896300           1054                    2549            0
#> 5     257.109   11897600           1469                    2526            0
#> 6     257.299   11899200           1795                    2507            0
#> 7     257.791   11901800           2250                    2636            0
#> 8     258.367   11903800           2643                    2798            0
#> 9     259.013   11904600           2987                    3019            0
#> 10    259.681   11909900           3256                    3280            0
#>    GVI.coverage alpha_i.gambiae_ss.
#> 1             0         0.000206347
#> 2             0         0.000206342
#> 3             0         0.000206349
#> 4             0         0.000206332
#> 5             0         0.000206344
#> 6             0         0.000206344
#> 7             0         0.000206347
#> 8             0         0.000206338
#> 9             0         0.000206330
#> 10            0         0.000206329
  1. The *_out.txt file
read.table("exp_test/outputs/exp_test_base_out.txt", header = FALSE, sep = "\t", nrows = 10)
#>    V1 V2 V3          V4
#> 1   1  1  0 6.07000e+02
#> 2   1  2  0 3.39300e+03
#> 3   1  1  3 1.37000e+02
#> 4   1  2  3 1.15500e+03
#> 5   1  1 14 3.35770e+04
#> 6   1  2 14 2.94240e+04
#> 7   1  1 15 7.11000e+02
#> 8   1  2 15 4.80000e+01
#> 9   1  0 35 2.19321e-01
#> 10  1  0 36 2.87284e-01

The *_out.txt file is usually the file you are interested in, as it contains the measurement data for each survey you defined in the input list. See here for more information.

Read results

OpenMalaria’s output files use numerical indiccators in the second column which can make the interpretation of the values difficult for new users. Thus OMU provides a helper function which reads an output file and assigns column names and translates certain values (e.g. the survey dates).

outfile <- readOutputFile("exp_test/outputs/exp_test_base_out.txt")
head(outfile)
#>    survey_date third_dimension measure value
#> 1:  1999-12-27               1   nHost   607
#> 2:  1999-12-27               2   nHost  3393
#> 3:  1999-12-27               1 nPatent   137
#> 4:  1999-12-27               2 nPatent  1155
#> 5:  1999-12-27               1 nUncomp 33577
#> 6:  1999-12-27               2 nUncomp 29424