Evaluate solution importance using incremental ranks

Calculate importance scores for planning units selected in a solution by calculating ranks via an incremental optimization process (based on Jung et al. 2021).

Usage

eval_rank_importance(x, solution, ...)

# S4 method for ConservationProblem,numeric
eval_rank_importance(x,
solution, ..., rescale, run_checks, force, by_zone, objective, extra_args,
n, budgets)

# S4 method for ConservationProblem,matrix
eval_rank_importance(x,
solution, ..., rescale, run_checks, force, by_zone, objective, extra_args,
n, budgets)

# S4 method for ConservationProblem,data.frame
eval_rank_importance(x,
solution, ..., rescale, run_checks, force, by_zone, objective, extra_args,
n, budgets)

# S4 method for ConservationProblem,Spatial
eval_rank_importance(x,
solution, ..., rescale, run_checks, force, by_zone, objective, extra_args,
n, budgets)

# S4 method for ConservationProblem,sf
eval_rank_importance(x, solution,
..., rescale, run_checks, force, by_zone, objective, extra_args, n, budgets)

# S4 method for ConservationProblem,Raster
eval_rank_importance(x,
solution, ..., rescale, run_checks, force, by_zone, objective, extra_args,
n, budgets)

# S4 method for ConservationProblem,SpatRaster
eval_rank_importance(x,
solution, ..., rescale, run_checks, force, by_zone, objective, extra_args,
n, budgets)

Arguments

x: problem() object.
solution: numeric, matrix, data.frame, terra::rast(), or sf::sf() object. The argument should be in the same format as the planning unit cost data in the argument to x. See the Solution format section for more information.
...: not used.
rescale: logical flag indicating if replacement cost values -- excepting infinite (Inf) and zero values -- should be rescaled to range between 0.01 and 1. Defaults to TRUE.
run_checks: logical flag indicating whether presolve checks should be run prior solving the problem. These checks are performed using the presolve_check() function. Defaults to TRUE. Skipping these checks may reduce run time for large problems.
force: logical flag indicating if an attempt should be made to solve the problem even if potential issues were detected during the presolve checks. Defaults to FALSE.
by_zone: logical indicating value. If TRUE, then the optimization process will increment budgets for each zone separately. If FALSE, then the optimization process will increment a single budget that is applied to all zones. Note that this parameter is only considered if n is specified, and does not affect processing if budgets is specified. Defaults to TRUE.
objective: character value with the name of the objective function that should be used for the incremental optimization process. This function must be budget limited (e.g., cannot be add_min_set_objective()). For example, "add_min_shortfall_objective" can be used to specify the minimum shortfall objective (per add_min_shortfall_objective()).. Defaults to NULL such that the same objective is used as specified in x. If using this default and x has the minimum set objective, then the minimum shortfall objective is used.
extra_args: list value with additional arguments for the objective function (excluding the budgets parameter). For example, this parameter can be used to supply phylogenetic data for the phylogenetic diversity objective function (i.e., when using objective = "add_max_phylo_div_objective"). Defaults to NULL such that no additional arguments are supplied.
n: integer number of ranks to evaluate. Note that either n or budgets (not both) must be specified. If n is specified, then by_zone is considered during processing for problems with multiple zones.
budgets: numeric vector with the budget thresholds for generating solutions at different steps in an iterative procedure. Note that either n or budgets (not both) must be specified.

Value

A numeric, matrix, data.frame, terra::rast(), or sf::sf() object containing the importance scores for each planning unit in the solution. Specifically, the returned object is in the same format as the planning unit data in the argument to x. The object also has the following attributes that provide information on the incremental optimization process.

budgets: numeric or matrix containing the budgets used for each increment in the incremental optimization process. If the problem (per x) has a single zone, then the budgets are a numeric vector, wherein values correspond to the budgets for each increment. Otherwise, if the problem (per x) has multiple zones, then the budgets are a matrix and their format depends on the by_zone parameter. If by_zone = FALSE, then the budgets are are matrix with a column for each zone and a row for each budget increment. Alternatively, if by_zone = TRUE, then the matrix has a single column and a row for each budget increment.
objective: numeric mathematical objective values for each solution generated during the incremental optimization process.
runtime: numeric total amount of time elapsed during the optimization (reported in seconds) of each solution generated throughout the incremental optimization process.
status: character status of the optimization process for each solution generated during the incremental optimization process. See solve() for further details.
gap: numeric optimality of each solution generated during the incremental optimization process. See solve() for further details.

Details

Importance scores are calculated using an incremental optimization process. Note that if a problem has complex constraints (i.e., constraints that do not involve locking in or locking out planning units), then the budgets parameter must be specified. This optimization process involves the following steps.

A set of budgets are defined. then the budgets are defined using the budgets. Otherwise, if an argument to the n parameter is supplied, then the budgets are defined as a set of n values with equal increments between them that sum to the total cost of solution. For example, if considering a problem with a single zone, a solution with a total cost of 400, and n = 4: then the budgets will be 100, 200, 300, and 400. If considering a multiple zone problem and by_zone = FALSE, then the budgets will based calculated based on the total cost of the solution across all zones. Otherwise if by_zone = TRUE, then the budgets are calculated and set based on the total cost of planning units allocated to each zone (separately) in the solution. Note that after running this function, you can see what budgets were used to calculate the ranks by accessing attributes from the result (see below for examples).
The problem (per x) is checked for potential issues. This step is performed to avoid issues during subsequent optimization steps. Note that this step can be skipped using run_checks = FALSE. Also, if issues are detected and you wish to proceed anyway, then useforce = TRUE ignore any detected issues.
The problem is modified for subsequent optimization. In particular, any planning units not selected in solution are locked out. This is important to ensure that all subsequent optimization procedures produce solutions that only contain planning units selected in the solution.
The problem is further modified for subsequent optimization. Specifically, its objective is overwritten using the objective defined for the rank calculations (per objective) with the smallest budget defined in the first step. Additionally, if an argument to the extra_args parameter is specified, this argument is also used when overwriting the objective.
The modified problem is solved to generate a solution. Depending on the budget and objective specified when modifying the problem, the newly generated solution will contain a subset of the planning units selected in the original solution. are assigned a rank. In particular, all selected planning units in the newly generated solution are assigned the same rank. This rank is based on the number of increments that have been previously completed. If no previous increments have been completed, then the rank is equal to the total number of budget increments. Otherwise, if previous increments have been completed, then the rank is equal to the total number of budget increments minus the number of completed increments. Note that if previous increments have been completed, then only planning units in the newly generated solution that have not been previously selected are assigned this rank. For example, if no previous increments have been completed and there are 5 budget increments (e.g. n = 5), then the planning units selected in the newly generated solution are assigned a rank of 5. Alternatively, if 3 previous increments have been completed, then the planning units would be assigned a rank of 2.
The problem is further modified for subsequent optimization. Specifically, the planning units selected in the newly generated solution are locked in. This is to ensure that all subsequent solutions will select these planning units and, in turn, build on this solution. Additionally, the newly generated solution is used to specify the starting solution for the subsequent optimization procedure to reduce processing time (note this is only done when using the CBC and Gurobi solvers).
Steps 4--7 are repeated for each of the remaining budget increments. As the increasingly greater budgets are used at higher increments, the modified problem will generate new solutions that become increasingly similar to the original solution. Planning units that are selected at lower budget increments are assigned greater ranks and considered more important, and those selected at higher budget increments are assigned lower ranks and considered less important. This is because if a planning unit is highly cost-effective for meeting the objective (per objective), then it is more likely to be selected earlier in the incremental optimization process.
The incremental optimization process has completed. If rescale = TRUE, then the ranks are linearly rescaled to range between 0.01 and 1. Otherwise, the ranks remain unchanged.
The rank values are output in the same format as the planning units in the problem (per x) (see the Solution Format section for details).

Solution format

Broadly speaking, the argument to solution must be in the same format as the planning unit data in the argument to x. Further details on the correct format are listed separately for each of the different planning unit data formats:

x has numeric planning units: The argument to solution must be a numeric vector with each element corresponding to a different planning unit. It should have the same number of planning units as those in the argument to x. Additionally, any planning units missing cost (NA) values should also have missing (NA) values in the argument to solution.
x has matrix planning units: The argument to solution must be a matrix vector with each row corresponding to a different planning unit, and each column correspond to a different management zone. It should have the same number of planning units and zones as those in the argument to x. Additionally, any planning units missing cost (NA) values for a particular zone should also have a missing (NA) values in the argument to solution.
x has terra::rast() planning units: The argument to solution be a terra::rast() object where different grid cells (pixels) correspond to different planning units and layers correspond to a different management zones. It should have the same dimensionality (rows, columns, layers), resolution, extent, and coordinate reference system as the planning units in the argument to x. Additionally, any planning units missing cost (NA) values for a particular zone should also have missing (NA) values in the argument to solution.
x has data.frame planning units: The argument to solution must be a data.frame with each column corresponding to a different zone, each row corresponding to a different planning unit, and cell values corresponding to the solution value. This means that if a data.frame object containing the solution also contains additional columns, then these columns will need to be subsetted prior to using this function (see below for example with sf::sf() data). Additionally, any planning units missing cost (NA) values for a particular zone should also have missing (NA) values in the argument to solution.
x has sf::sf() planning units: The argument to solution must be a sf::sf() object with each column corresponding to a different zone, each row corresponding to a different planning unit, and cell values corresponding to the solution value. This means that if the sf::sf() object containing the solution also contains additional columns, then these columns will need to be subsetted prior to using this function (see below for example). Additionally, the argument to solution must also have the same coordinate reference system as the planning unit data. Furthermore, any planning units missing cost (NA) values for a particular zone should also have missing (NA) values in the argument to solution.

References

Jung M, Arnell A, de Lamo X, García-Rangel S, Lewis M, Mark J, Merow C, Miles L, Ondo I, Pironon S, Ravilious C, Rivers M, Schepaschenko D, Tallowin O, van Soesbergen A, Govaerts R, Boyle BL, Enquist BJ, Feng X, Gallagher R, Maitner B, Meiri S, Mulligan M, Ofer G, Roll U, Hanson JO, Jetz W, Di Marco M, McGowan J, Rinnan DS, Sachs JD, Lesiv M, Adams VM, Andrew SC, Burger JR, Hannah L, Marquet PA, McCarthy JK, Morueta-Holme N, Newman EA, Park DS, Roehrdanz PR, Svenning J-C, Violle C, Wieringa JJ, Wynne G, Fritz S, Strassburg BBN, Obersteiner M, Kapos V, Burgess N, Schmidt- Traub G, Visconti P (2021) Areas of global importance for conserving terrestrial biodiversity, carbon and water. Nature Ecology and Evolution, 5: 1499--1509.

Examples

# \dontrun{
# seed seed for reproducibility
set.seed(600)

# load data
sim_pu_raster <- get_sim_pu_raster()
sim_pu_polygons <- get_sim_pu_polygons()
sim_features <- get_sim_features()
sim_zones_pu_raster <- get_sim_zones_pu_raster()
sim_zones_features <- get_sim_zones_features()

# create minimal problem with binary decisions
p1 <-
  problem(sim_pu_raster, sim_features) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.1) %>%
  add_binary_decisions() %>%
  add_default_solver(gap = 0, verbose = FALSE)

# solve problem
s1 <- solve(p1)

# print solution
print(s1)
#> class       : SpatRaster 
#> dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
#> resolution  : 0.1, 0.1  (x, y)
#> extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#> coord. ref. : Undefined Cartesian SRS 
#> source(s)   : memory
#> varname     : sim_pu_raster 
#> name        : layer 
#> min value   :     0 
#> max value   :     1 

# plot solution
plot(s1, main = "solution", axes = FALSE)


# calculate importance scores using ranks based on 10 budgets
# N.B. since the objective for calculating ranks is not explicitly
# defined and the problem has a minimum set objective, the
# ranks are calculated using the minimum shortfall objective by default
rs1 <- eval_rank_importance(p1, s1, n = 10)

# print importance scores
print(rs1)
#> class       : SpatRaster 
#> dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
#> resolution  : 0.1, 0.1  (x, y)
#> extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#> coord. ref. : Undefined Cartesian SRS 
#> source(s)   : memory
#> varname     : sim_pu_raster 
#> name        : rs 
#> min value   :  0 
#> max value   :  1 

# plot importance scores
plot(rs1, main = "rank importance (10, min shortfall obj", axes = FALSE)


# display optimization information from the attributes
## status
print(attr(rs1, "status"))
#>  [1] "OPTIMAL" "OPTIMAL" "OPTIMAL" "OPTIMAL" "OPTIMAL" "OPTIMAL" "OPTIMAL"
#>  [8] "OPTIMAL" "OPTIMAL" "OPTIMAL"
## optimality gap
print(attr(rs1, "gap"))
#>  [1] 0 0 0 0 0 0 0 0 0 0
## run time
print(attr(rs1, "runtime"))
#>  [1] 0.005 0.003 0.003 0.002 0.003 0.002 0.003 0.003 0.003 0.002
## objective value
print(attr(rs1, "objective"))
#>  [1] 4.4831422 3.9636924 3.4483566 2.9239906 2.4229403 1.8946389 1.3955220
#>  [8] 0.8733925 0.3850706 0.0000000

# plot relationship between objective values and rank
plot(
  y = attr(rs1, "objective"),
  x = seq_along(attr(rs1, "objective")),
  ylab = "objective value", xlab = "rank",
  main = "relationship between objective values and rank"
)


# calculate importance scores using the maximum utility objective and
# based on 10 different budgets
rs2 <- eval_rank_importance(
  p1, s1, n = 10, objective = "add_max_utility_objective"
)
#> Warning: Targets specified for the problem will be ignored.
#> ℹ If the targets are important, use a different objective.

# print importance scores
print(rs2)
#> class       : SpatRaster 
#> dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
#> resolution  : 0.1, 0.1  (x, y)
#> extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#> coord. ref. : Undefined Cartesian SRS 
#> source(s)   : memory
#> varname     : sim_pu_raster 
#> name        : rs 
#> min value   :  0 
#> max value   :  1 

# plot importance scores
plot(rs2, main = "rank importance (10, max utility obj)", axes = FALSE)


# calculate importance scores using ranks based on 5 manually specified
# budgets

# calculate 5 ranks using equal intervals
# N.B. we use length.out = 6 because we want 5 budgets > 0
budgets <- seq(0, eval_cost_summary(p1, s1)$cost[[1]], length.out = 6)[-1]

# calculate importance using manually specified budgets
# N.B. since the objective for calculating ranks is not explicitly
# defined and the problem has a minimum set objective, the
# ranks are calculated using the minimum shortfall objective by default
rs3 <- eval_rank_importance(p1, s1, budgets = budgets)

# print importance scores
print(rs3)
#> class       : SpatRaster 
#> dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
#> resolution  : 0.1, 0.1  (x, y)
#> extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#> coord. ref. : Undefined Cartesian SRS 
#> source(s)   : memory
#> varname     : sim_pu_raster 
#> name        : rs 
#> min value   :  0 
#> max value   :  1 

# plot importance scores
plot(rs3, main = "rank importance (manual)", axes = FALSE)


# build multi-zone conservation problem with raster data
p4 <-
  problem(sim_zones_pu_raster, sim_zones_features) %>%
  add_min_set_objective() %>%
  add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>%
  add_binary_decisions() %>%
  add_default_solver(gap = 0, verbose = FALSE)

# solve the problem
s4 <- solve(p4)
names(s4) <- paste0("zone ", seq_len(terra::nlyr(sim_zones_pu_raster)))

# print solution
print(s4)
#> class       : SpatRaster 
#> dimensions  : 10, 10, 3  (nrow, ncol, nlyr)
#> resolution  : 0.1, 0.1  (x, y)
#> extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#> coord. ref. : Undefined Cartesian SRS 
#> source(s)   : memory
#> varnames    : sim_zones_pu_raster 
#>               sim_zones_pu_raster 
#>               sim_zones_pu_raster 
#> names       : zone 1, zone 2, zone 3 
#> min values  :      0,      0,      0 
#> max values  :      1,      1,      1 

# plot solution
# each panel corresponds to a different zone, and data show the
# status of each planning unit in a given zone
plot(s4, axes = FALSE)


# calculate importance scores
rs4 <- eval_rank_importance(p4, s4, n = 5)
names(rs4) <- paste0("zone ", seq_len(terra::nlyr(sim_zones_pu_raster)))

# plot importance
# each panel corresponds to a different zone, and data show the
# importance of each planning unit in a given zone
plot(rs4, axes = FALSE)

# }