Sensitivity analysis of crop models such as APSIM is a valuable task which allows us to develop an intuition of the model response to changing a variety of inputs. These inputs can generally be included in the following categories:
The weather component can depend on different sources of data or manipulations to simulate different conditions and stresses. The soil profiles can be from different locations or specifically manipulated to understand the impact of particular soil properties. Crop properties or cultivar specific parameters can also be investigated as a source of model behavior. Finally, examples of management usually include planting date, fertilizer application, tillage or residue management and irrigation. All of these components can be sources of uncertainty and can be evaluated in a sensitivity analysis.
When multiple parameters are evaluated performing model runs for all combinations can become prohibitive. For this reason, there are some combinations which can be more efficient and provide nearly the same amount of information. The suggestion here is to build complexity in the sensitivity analysis task gradually. For example, first choose one parameter and select 5-10 reasonable values, run the model and evaluate the model output. This is a one parameter at a time approach and it can be of great value in developing an intuition for the model behavior.
For a general introduction to sensitivity analysis see Saltelli et al. 2008. Global sensitivity analysis. The Primer.
Here they define sensitivity analysis as:
The study of how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input
A preliminary step in sensitivity analysis often involves uncertainty analysis and, as an example, it can involve computing the output from a model based on various samples from the relevant inputs. Usually, this involves identifying parameters, choosing appropriate distributions which we sample from and creating candidate values which are then used in the model.
The steps needed to perform a uncertainty analysis can involve simple use of for loops.
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.
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).
parm.vector.index = some paths point to a vector of parameters, but sometimes only one element of that vector needs to be evaluated. In this case provide the index for which element of the vector needs to be evaluated. If some of the parameters need this index while others do not, include a zero for those parameter vectors which are not indexed.
replacement = logical indicating whether each parameter is in the replacement component or not.
grid = data frame with combinations of parameters to be evaluated.
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.
… = additional arguments (none used at the moment).
The example below shows the following steps.
Here the grid is created using expand.grid and it has 9 rows. This means that this while 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)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.
## 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.0A 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))## Warning in morris(factors = c("Fertiliser", "Population"), r = 3, design =
## list(type = "oat", : keeping 2 repetitions out of 3## Fertiliser Population
## [1,] 100 5.0
## [2,] 200 5.0
## [3,] 200 152.5
## [4,] 200 5.0
## [5,] 200 152.5
## [6,] 100 152.5The 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.