In JGCRI/gcamland: GCAM Land Allocation Module

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

knitr::opts_chunk$set(echo = TRUE)
library(gcamland)
library(dplyr)
library(tidyr)

Setup params

Pick the number of runs to do. We select only two for this tutorial to speed things up.

Nruns <- 2

Do some runs

By default, no analysis is run. Only the ensemble of runs is generated, and analysis can be done after the fact.

An example of the call to run the model, as well as the progress print statements included.

set.seed(1)

scenObjects <- run_ensemble(N=Nruns)

Run the objective function analysis

• Right now, the run_objective function evaluates user-specified objective function for each land type separately. Certainly can give an argument with fewer land types.
• idea is that it's all so fast and not that big, just go ahead and calculate for each land type. Then when we call minimize_objective, that function also takes a list of land types of interest, averages the objective function across those land types, and then picks the parameter set that gives you the smallest average objective function value (ie minimized error) across those land types.

• Also, currently, from get_historical_land_data in utility-functions-analysis.R (previously in bayesian.R), this only supports comparison to FAO GCAM commodities - mostly crop data and not Forest or anything. A note to change this is included in the documentation of minimize_objective.

• Note that land cover types X region combinations that have a standard deviation of 0 in the observational data, warnings will be thrown and the metrics nrms, centerednrms, and kge will return NA values. While we have not seen this in observational data in the United States, we have not yet checked every commodity in every region.

• For land cover types X region combinations that have a standard deviation of 0 in the modeled data, warnings will be thrown and the metric kge will return NA values. In the United States, this means PalmFruit.

An example of the function call to run the objective function analysis.

scenObjectsEvaluated <- run_objective(scenObjects)

And different views of the results of that analysis. Note that the list item $mObjFunEval for each scenario is where this analysis is contained:

str(scenObjectsEvaluated)

names(scenObjectsEvaluated[[1]])

knitr::kable(head(scenObjectsEvaluated[[1]]$mObjFunEval))

knitr::kable(tail(scenObjectsEvaluated[[1]]$mObjFunEval))

knitr::kable(head(scenObjectsEvaluated[[1]]$mObjFunEval %>%
                    filter(land.type == 'PalmFruit', 
                           objfun == 'kge')))

Organize these evaluated results

Using the grand_table_objective function results of the objective function evaluations are reorganized into a grand table of parameters. The table produced includes the model parameters and all user-specified objective function evaluations, for all of the models in the input list.

An example of the call to this function and a quick look at the structure of the results.

grand_table_objective(aScenarioList = scenObjectsEvaluated) %>%
  # these lines of code are because the `left_join` in 
  # `add_parameter_data` (from utility-functions-analysis.R) is only
  # a join on scenario, rather than scenario and region.
  #
  # TODO:
  # Unclear if should update that join or if we somehow want to keep it
  # separate; ideally will figure out more when do other regions.
  # For now, leaving it alone so that the `grand_table` works as it
  # did before. 
  # (relevant because both `grand_table` and `grand_table_objective`
  # call `add_parameter_data`)
  filter(region.x == region.y) %>%
  rename(region = region.x) %>%
  select(-region.y) ->
  GTobjective

knitr::kable(head(GTobjective))

And a view of a reshaped, 'prettier' (that is, wide) version of these results. However, the long version GTobjective is what the next step in analysis is expecting.

GTobjective %>%
  spread(objfun, objfunval) ->
  GTobjective_pretty

knitr::kable(head(GTobjective_pretty))

Select the parameter set from runs that minimize specified obj fcn across specified land types

minimize_objective is a minimizer function. Itcans the runs in a given sample of model runs and selects the parameter set from those runs that minimizes the mean across user-specified landtypes of a user-specified objective function, objfun_to_min.

Exceptions

For bias, the quantity being minimized is actually the average over specified land types of abs(bias), as a bias of 0 is perfect performance, and a bias of 0.25 vs -0.25 are equally 'good'.

For kge, the quantity being minimized is actually the average over specified land types of 1 -kge. This is chosen because of the following:

1. Values of kge can range anywhere from -infinity to +1, with a value of +1 being perfect model performance, therefore kge is not directly useful as a quantity to be minimized.

2. Values of 1-kge range from 0 to +infinity, with 0 corresponding to perfect model performance, as well as being the minimum value achieved by this measure. Additionally, because all values of this quantity are positive, no cancelation of errors may occur when calculating the mean across specified land types.

Further details are given in the documentation of the minimize_objective function.

default arguments

Default arguments of minimize_objective currently is focused on RMS (objfun_to_min = 'rms') for all GCAM commodities landtypes = c( "Corn", "FiberCrop", "MiscCrop", "OilCrop", "OtherGrain", "PalmFruit", "Rice", "Root_Tuber", "SugarCrop", "Wheat") that we run and have FAO comparison data for.

minimize_objective is taking the average of RMS across all of these land types, and selecting the parameter set that minimizes that quantity.

Currently no support for mixing and matching different objective functions for different crops. Seems doable but not implemented yet (eg maybe you want parameters that minimize both Corn's bias and OtherGrain's KGE.)

An example call to the minimize_objective function assuming default arguments, and a printout of the results:

knitr::kable(minimize_objective(GTobjective))

Corn, OilCrop bias example

Select parameters that minimize the average (absolute value of) bias in Corn and OilCrop.

The current implementation does not allow for cancellation of errors: the average across specified land types of abs(bias) is what is actually being minimized in minimize_objective. So, for example, a positive bias in Corn and a negative bias in OilCrop will not cancel out to give an artifically small mean error.

knitr::kable(minimize_objective(GTobjective, objfun_to_min = 'bias', landtypes = c('Corn', 'OilCrop')))

NOTE

The error reported in the column landTypeMeanObjFunVal is what you think: it's the average error only across the specified land types. If you want to see what the error on Wheat is for this example, the results from grand_table_objective (ie GTobjective) must be filtered and examined.

GTobjective %>%
   filter(land.type == 'Wheat')

OtherGrain rms example

knitr::kable(minimize_objective(GTobjective, objfun_to_min = 'rms', landtypes = c('OtherGrain')))

JGCRI/gcamland documentation built on Oct. 6, 2020, 5:30 p.m.