Uncertainty Analysis using for loops

The steps needed to perform a uncertainty analysis can involve simple use of for loops.

1. Identify a parameter of interest
2. Sample (or select) a new value for the parameter of interest
3. Edit the APSIM file with the new value
4. Run the model
5. Collect the results and repeat as needed

A simple example

## Create a temporary directory and copy the Wheat.apsimx file
tmp.dir <- tempdir()
extd.dir <- system.file("extdata", package = "apsimx")
file.copy(file.path(extd.dir, "Wheat.apsimx"), tmp.dir)
## Identify a parameter of interest
## In this case we want to know the impact of varying the fertilizer amount
pp <- inspect_apsimx("Wheat.apsimx", src.dir = tmp.dir, 
                      node = "Manager", parm = list("SowingFertiliser", 1))
## For simplicity, we create a vector of fertilizer amounts (instead of sampling)
ferts <- seq(5, 200, length.out = 7)
col.res <- NULL
for(i in seq_along(ferts)){
  
  edit_apsimx("Wheat.apsimx", src.dir = tmp.dir,
              node = "Other",
              parm.path = pp, parm = "Amount",
              value = ferts[i])
  
  sim <- apsimx("Wheat.apsimx", src.dir = tmp.dir)
  col.res <- rbind(col.res, data.frame(fertilizer.amount = ferts[i], 
                                       wheat.aboveground.wt = mean(sim$Wheat.AboveGround.Wt, na.rm = TRUE)))
}

The code here is shown for illustration so that the basic strategy can be replicated. The results collected in object col.res can be further processed to better understand the impact of the chosen parameter (input) on the output.

There is a function in the package which simplifies this task. The one below is for APSIM Next Gen and there is a similar one for APSIM Classic.

sens_apsimx

• file = .apsimx file

• src.dir = source directory where the APSIM file is located.

• parm.paths = character vector with the full path to the parameters to be evaluated. This can be obtained by using inspect_apsimx or inspect_apsim_json. The length of the parameter vector should be equal to the number of parameters being evaluated (sometimes less is more).

• convert = whether to convert strings to numeric

• replacement = logical indicating whether each parameter is in the replacement component or not.

• grid = data frame with combinations of parameters to be evaluated.

• soil.profiles = soil profiles, which can be passed as a list

• summary = optional summary function to apply if the APSIM output returns multiple rows. The mean is used by default.

• root = indicate the root simulation in case it is not the default one.

• verbose = printing progress and other details

• cores = number of cores to use. The package ‘future’ is currently used.

• save = whether to save intermediate (and final) results.

• … = additional arguments (none used at the moment).

Example using sens_apsimx

The example below shows the following steps.

1. Copy the example ‘Wheat’ file to at temporary directory
2. Identify parameters of interest (using inspect_apsimx)
3. Create a simple grid for evaluating the model

Here the grid is created using expand.grid and it has 9 rows. This means that this will require 9 simulation runs of the APSIM model. This takes a couple of minutes to run. The resulting object is of class sens_apsim and it can be investigated using the summary function (for now).

tmp.dir <- tempdir()
extd.dir <- system.file("extdata", package = "apsimx")
file.copy(file.path(extd.dir, "Wheat.apsimx"), tmp.dir)
## Identify a parameter of interest
## In this case we want to know the impact of varying the fertilizer amount
## and the plant population
pp1 <- inspect_apsimx("Wheat.apsimx", src.dir = tmp.dir, 
                      node = "Manager", parm = list("SowingFertiliser", 1))
pp1 <- paste0(pp1, ".Amount")

pp2 <- inspect_apsimx("Wheat.apsimx", src.dir = tmp.dir, 
                      node = "Manager", parm = list("SowingRule1", 9))
pp2 <- paste0(pp2, ".Population")

## The names in the grid should (partially) match the parameter path names
grd <- expand.grid(Fertiliser = c(50, 100, 150), Population = c(100, 200, 300))
 
## This takes 2-3 minutes
sns <- sens_apsimx("Wheat.apsimx", src.dir = tmp.dir,
                    parm.paths = c(pp1, pp2),
                    grid = grd)

Creating grids

There are many ways of creating grids. The example above using function expand.gird returned all possible combinations of the parameter values. This might not be desirable or feasible and other methods are available. For this task, the package sensitivity is especially useful.

library(sensitivity)

## Registered S3 method overwritten by 'sensitivity':
##   method    from 
##   print.src dplyr

## Simple example: see documentation for other options
X.grid <- parameterSets(par.ranges = list(Fertiliser = c(1, 300), 
                                          Population = c(1, 300)),
                        samples = c(3,3), method = "grid")
X.grid

##       Fertiliser Population
##  [1,]        1.0        1.0
##  [2,]      150.5        1.0
##  [3,]      300.0        1.0
##  [4,]        1.0      150.5
##  [5,]      150.5      150.5
##  [6,]      300.0      150.5
##  [7,]        1.0      300.0
##  [8,]      150.5      300.0
##  [9,]      300.0      300.0

Second option: Morris

A second option is to use the function morris.

X.mrrs <- morris(factors = c("Fertiliser", "Population"),
                   r = 3, design = list(type = "oat", levels = 3, grid.jump = 1),
                   binf = c(0, 5), bsup = c(200, 300))
X.mrrs$X

##       Fertiliser Population
##  [1,]        100      152.5
##  [2,]        100        5.0
##  [3,]        200        5.0
##  [4,]        100        5.0
##  [5,]        100      152.5
##  [6,]          0      152.5
##  [7,]          0      152.5
##  [8,]          0        5.0
##  [9,]        100        5.0

The sensitivity package allows for decoupling the simulation runs and the analysis of the results. For example:

## This takes 2-3 minutes
sns2 <- sens_apsimx("Wheat.apsimx", src.dir = tmp.dir,
                     parm.paths = c(pp1, pp2),
                     grid = X.mrrs$X)
## These are the sensitivity results for AboveGround.Wt only
sns.res.ag <- tell(X.mrrs, sns2$grid.sims$Wheat.AboveGround.Wt)
sns.res.ag
## Call: morris(factors = c("Fertiliser", "Population"), 
## r = 3, design = list(type = "oat", levels = 3, grid.jump = 1), 
## binf = c(0, 5), bsup = c(200, 300))
## 
## Model runs: 9 
##                 mu  mu.star    sigma
## Fertiliser 916.6674 916.6674 662.8798
## Population 448.4530 448.4530 640.6807
plot(sns.res.ag)

I’m not running the code in this vignette to save time, but hopefully this is sufficient information to illustrate performing sensitivity analysis using package apsimx and sensitivity.