Decision-makers often wish to have a quantitative basis for their decisions. However, there is often no ‘hard data’ for many important variables, which can paralyze decision-making processes or lead decision-makers to conclude that large research efforts are needed before a decision can be made. That is, many variables decision makers must consider cannot be precisely quantified, at least not without unreasonable effort. The major objective of (prescriptive) decision analysis is to support decision-making processes faced with this problem (E. Luedeling and Shepherd 2016). Decision analysis can make forecasts of decision outcomes without precise numbers, as long as probability distributions describing the possible values for all variables can be estimated.
The decisionSupport package (Eike Luedeling and Göhring 2017) implements this as a Monte Carlo simulation, which generates a large number of plausible system outcomes, based on random numbers for each input variable that are drawn from user-specified probability distributions. It also conducts a sensitivity analysis (based on Projection-to-Latent-Structures regression) to highlight uncertain variables that have large impacts on model outputs (Eike Luedeling and Gassner 2012; Wold, Sjostrom, and Eriksson 2001). This approach is useful for determining whether a clearly preferable course of action can be delineated based on the present state of knowledge without the need for further information. If the distribution of predicted system outcomes doesn’t imply a clearly preferable decision option, variables identified as carrying decision-relevant uncertainty can then be targeted by decision-supporting research (Eike Luedeling et al. 2015). This approach is explained in more detail below.
The decisionSupport() function in the R package decisionSupport can be applied to conduct decision analysis. The function requires two inputs:
These two inputs are provided as arguments to the decisionSupport() function, which conducts a Monte Carlo analysis with repeated model runs based on probability distributions for all uncertain variables. The data table and model are customized to fit the particulars of a specific decision.
The ecology of the conifer forests of the American west is largely shaped by wildfires. Many ecologists recommend regular and low-intensity burns to manage the build-up of combustible understory materials (shrubs, fallen branches). However, not all municipalities or regions implement these practices. Failure to follow the recommended controlled burning regime may lead to the build-up of fire stock and increase the severity of wildfires. Such wildfires have destroyed many conifer forest ecosystems of the Western United States in the recent past.
Here, we provide a simple example (in annotated R code) of how the decisionSupport package can be used to inform a decision process. The example provided simulates the decision of forest managers to use controlled fires in conifer forests vs. running the risk of severe fire.
To support the model building process we design an input table wildfire_input_table.csv containing some of the basic values for the analysis. This table contains all the variables used in the model. Their distributions are described by a 90% confidence interval, which is specified by lower (5% quantile) and upper (95% quantile) bounds, as well as the shape of the distribution. This example uses four different distributions:
const – a constant valuenorm – a normal distributiontnorm_0_1 – a truncated normal distribution that can only have values between 0 and 1 (useful for probabilities; note that 0 and 1, as well as numbers outside this interval are not permitted as inputs)posnorm – a normal distribution truncated at 0 (only positive values allowed)Here’s the input table for the wildfire example:
setwd(getwd()) # this should be replaced by the folder where the wildfire_input_table.csv
#file it stored
format(read.csv("wildfire_input_table.csv")[,1:5],scientific=FALSE)## variable lower median upper distribution
## 1 n_years 20.0 NA 20.00 const
## 2 area_size 1.0 NA 1.00 const
## 3 initial_biomass 20.0 NA 100.00 posnorm
## 4 biomass_after_burning 10.0 NA 20.00 posnorm
## 5 biomass_after_severe_fire 2.0 NA 5.00 posnorm
## 6 biomass_after_mild_fire 10.0 NA 20.00 posnorm
## 7 bio_accu_rate_contr_burn 10.0 NA 40.00 posnorm
## 8 bio_accu_rate_no_contr_burn 40.0 NA 100.00 posnorm
## 9 severity_threshold 50.0 NA 120.00 posnorm
## 10 discount_rate 2.0 NA 6.00 posnorm
## 11 controlled_burning_frequency 2.0 NA 2.00 const
## 12 fire_risk 0.1 NA 0.25 tnorm_0_1
## 13 cost_of_controlled_burning 300.0 NA 1000.00 posnorm
## 14 env_imp_severe_fire 1000.0 NA 6000.00 posnorm
## 15 env_imp_mild_fire -500.0 NA 800.00 norm
## 16 fire_fighting_cost_mild_fire 100.0 NA 500.00 posnorm
## 17 fire_fighting_cost_severe_fire 1000.0 NA 10000.00 posnormFor a full list of possible distributions, type ?random.estimate1d in R. When specifying confidence intervals for truncated distributions, note that approximately 5% of the random values should ‘fit’ within the truncation interval on either side. If there isn’t enough space, the function will generate a warning (usually it will still work, but the inputs may not look like you intended them to).
Default distributions are provided for all variables, but feel free to make adjustments by editing the .csv file in a spreadsheet program.
To set up the analysis we first define a time line in the input table that contains information about whether there is a wildfire or a controlled burn in each year of our planning horizon. We specify this using the variable n_years (we use 20 years in this example).
We define a probability fire_risk (which can range from 0 to 1) that a wildfire will break out in a given year, and use the chance_event() function to simulate a time series of wildfire occurrence:
Fire <- chance_event(fire_risk, value_if = 1, n = n_years)This line of R code produces a series of n_years (here: 20) numbers (called a vector in R), which have a fire_risk chance of being 1 (specified by the value_if argument) and are 0 otherwise. The 1s indicate years that feature a wildfire, 0s mark years without a fire.
We also simulate whether controlled burns occur in a given year. If we assume that fire managers follow a strict schedule, this can be computed as
Controlled <- as.numeric(1:n_years %% controlled_burning_frequency==0)In the above code controlled_burning_frequency is from the input table, which is simply the number of years between two controlled burns (should be an integer). The symbol %% calls a modulo operation, which returns the remainder after dividing each of the year numbers (for 1:n_years) by controlled_burning_frequency. If the year number is a multiple of the controlled_burning_frequency the code returns a TRUE evaluation. If the year is not a multiple, it returns a FALSE. These values are then converted into 0s and 1s by using the as.numeric() function.
The code below combines these vectors into a data.frame called sim.fire that contains all information on our fire scenario.
sim.fire <- data.frame(Year=1:n_years, Fire, Controlled)
sim.fire## Year Fire Controlled
## 1 1 0 0
## 2 2 1 1
## 3 3 0 0
## 4 4 0 1
## 5 5 1 0
## 6 6 0 1
## 7 7 0 0
## 8 8 0 1
## 9 9 0 0
## 10 10 0 1
## 11 11 1 0
## 12 12 1 1
## 13 13 0 0
## 14 14 0 1
## 15 15 0 0
## 16 16 1 1
## 17 17 1 0
## 18 18 0 1
## 19 19 0 0
## 20 20 0 1We want to distinguish between mild and severe fires, because these differ strongly in the impact on the ecosystem and in the costs incurred in fighting the fire. To determine whether a fire is mild or severe, we simulate the amount of combustible biomass that has accumulated in a given year. When this biomass exceeds the threshold severity_threshold (from the input table) at the time of a wildfire outbreak, we classify the fire as severe, otherwise as mild.
To do this, we first have to simulate the amount of combustible biomass. We start by setting up two data.frames to store the results of this simulation, one for the no-burn scenario (no.control), and one for the scenario with controlled burns (control):
no.control <- data.frame(Year = integer(),
Biomass = double(),
FireType = character(),
stringsAsFactors = F)
control <- data.frame(Year = integer(),
Biomass = double(),
FireType = character(),
stringsAsFactors = F)Each of the two data.frames has three columns for the simulation year (Year), the combustible biomass (Biomass) and the fire type (FireType). For year 1, we set the biomass for both scenarios to an initial level described by initial_biomass from the input table. We then start a loop that cycles through all years (with y used as a counter to describe which year we’re in) and computes that year’s biomass, according to the fire scenario that is encountered in each year. We distinguish between four different fire settings (using the switch function): 1) When neither fire nor controlled burns occur, biomass increases by a certain percentage (bio_accu_rate_no_contr_burn). 2) When a controlled burn occurs, biomass is reduced to biomass_after_burning, but then increases by a certain percentage (bio_accu_rate_contr_burn). 3) When a mild fire occurs, the biomass is reduced to biomass_after_mild_fire. 4) When a severe fire occurs, the biomass declines to biomass_after_severe_fire.
All variables in steps 1-4 are from the input table.
We can then determine the fire severity by comparing the combustible biomass level with the fire severity. For this, we use several ifelse calls that first check if there was a controlled burn, then (if not) whether there was a fire, and finally whether biomass exceeded the threshold for a severe fire. Depending on the outcome of these tests, the FireType is set to “Controlled,” “None,” “Severe” or “Mild.”
This is simulated with the following code for the control scenario:
#start of year loop
for (y in 1:n_years) {
# record year
control[y, "Year"] <- y
# calculate biomass
if (y == 1) {control[y, "Biomass"] <- initial_biomass} else {
switch(EXPR = control[y - 1, "FireType"],
"Controlled" = {control[y, "Biomass"] <- biomass_after_burning *
(1 + bio_accu_rate_contr_burn/100)},
"None" = {control[y, "Biomass"] <- control[y - 1, "Biomass"] *
(1 + bio_accu_rate_no_contr_burn/100)},
"Severe" = {control[y, "Biomass"] <- biomass_after_severe_fire},
"Mild" = {control[y, "Biomass"] <- biomass_after_mild_fire})}
# determine wild fire type
ifelse(test = sim.fire[y,"Controlled"] == 1,
yes = control[y,"FireType"] <- "Controlled",
no = ifelse(test = sim.fire[y,"Fire"] == 0,
yes = control[y,"FireType"] <- "None",
no = ifelse(test = sim.fire[y,"Fire"] == 1 && control[y,"Biomass"] > severity_threshold,
yes = control[y,"FireType"] <- "Severe",
no = control[y,"FireType"] <- "Mild")))
} #end of year loop controlled burning scenarioFor the no.control scenario, the code looks similar, except that there are no consequences of a Controlled fire, because such fires don’t occur here:
#start of year loop
for (y in 1:n_years) {
# record year
no.control[y, "Year"] <- y
# calculate biomass
if (y == 1) {no.control[y, "Biomass"] <- initial_biomass} else {
switch(EXPR = no.control[y - 1, "FireType"],
"None" = {no.control[y, "Biomass"] <- no.control[y - 1, "Biomass"] *
(1 + bio_accu_rate_no_contr_burn/100)},
"Severe" = {no.control[y, "Biomass"] <- biomass_after_severe_fire},
"Mild" = {no.control[y, "Biomass"] <- biomass_after_mild_fire})
}
# determine fire type
ifelse(test = sim.fire[y,"Fire"] == 0,
yes = no.control[y,"FireType"] <- "None",
no = ifelse(
test = no.control[y,"Biomass"] > severity_threshold,
yes = no.control[y,"FireType"] <- "Severe",
no = no.control[y,"FireType"] <- "Mild"))
} #end of year loop 'no burn' scenarioAt this point, the fire time series is complete, with all information stored in the data.frames control and no.control:
control## Year Biomass FireType
## 1 1 43.16919 None
## 2 2 63.65350 Controlled
## 3 3 15.67209 None
## 4 4 23.10869 Controlled
## 5 5 15.67209 Mild
## 6 6 20.92580 Controlled
## 7 7 15.67209 None
## 8 8 23.10869 Controlled
## 9 9 15.67209 None
## 10 10 23.10869 Controlled
## 11 11 15.67209 Mild
## 12 12 20.92580 Controlled
## 13 13 15.67209 None
## 14 14 23.10869 Controlled
## 15 15 15.67209 None
## 16 16 23.10869 Controlled
## 17 17 15.67209 Mild
## 18 18 20.92580 Controlled
## 19 19 15.67209 None
## 20 20 23.10869 Controlledno.control## Year Biomass FireType
## 1 1 43.169186 None
## 2 2 63.653500 Severe
## 3 3 3.310227 None
## 4 4 4.880971 None
## 5 5 7.197053 Mild
## 6 6 20.925804 None
## 7 7 30.855357 None
## 8 8 45.496607 None
## 9 9 67.085312 None
## 10 10 98.918125 None
## 11 11 145.856002 Severe
## 12 12 3.310227 Mild
## 13 13 20.925804 None
## 14 14 30.855357 None
## 15 15 45.496607 None
## 16 16 67.085312 Severe
## 17 17 3.310227 Mild
## 18 18 20.925804 None
## 19 19 30.855357 None
## 20 20 45.496607 NoneIn the above simulation, the controlled_burning_frequency was 2 and the severity_threshold 53.
With this information we can move on to computing the cost implications for each scenario. We only consider three types of costs: the cost of controlled burning (cost_of_controlled_burning), the cost of fighting fires (fire_fighting_cost), and the value of the environmental damage caused by a fire (env_imp). We distinguish between severe and mild fires, so we get separate variables for the cost of firefighting (fire_fighting_cost_mild_fire and fire_fighting_cost_severe_fire) and the environmental impacts (env_imp_mild_fire and env_imp_severe_fire). For each year and each scenario, the applicable costs are selected (depending on the fire type) and added up to produce the total cost for each option. The result is then multiplied by the area under consideration (area_size) to produce the total cost.
# calculate costs of no.control treatment -------------------------
ff.cost.no.control <-
(no.control$FireType == "Severe") * fire_fighting_cost_severe_fire +
(no.control$FireType == "Mild") * fire_fighting_cost_mild_fire
env.cost.no.control <-
(no.control$FireType == "Severe") * env_imp_severe_fire +
(no.control$FireType == "Mild") * env_imp_mild_fire
bottomline_no_burn <- - ( ff.cost.no.control + env.cost.no.control) * area_size
# calculate costs of control treatment
treatment.cost.control <-
(control$FireType == "Controlled") * cost_of_controlled_burning
ff.cost.control <-
(control$FireType == "Severe") * fire_fighting_cost_severe_fire +
(control$FireType == "Mild") * fire_fighting_cost_mild_fire
env.cost.control <-
(control$FireType == "Severe") * env_imp_severe_fire +
(control$FireType == "Mild") * env_imp_mild_fire
bottomline_contr_burn <-
- (treatment.cost.control + ff.cost.control + env.cost.control) * area_sizeNow we have two vectors that contain the annual costs of the two scenarios (bottomline_contr_burn and bottomline_no_burn). The net benefits of the controlled burning strategy (benefit_of_controlled_burning) are then simply calculated as the difference between control and no control:
benefit_of_controlled_burning <- bottomline_contr_burn - bottomline_no_burnThis produces, for example:
bottomline_no_burn## [1] 0.0000 -7279.2190 0.0000 0.0000 51.9453 0.0000
## [7] 0.0000 0.0000 0.0000 0.0000 -7279.2190 51.9453
## [13] 0.0000 0.0000 0.0000 -7279.2190 51.9453 0.0000
## [19] 0.0000 0.0000bottomline_contr_burn## [1] 0.0000 -786.5574 0.0000 -786.5574 51.9453 -786.5574 0.0000
## [8] -786.5574 0.0000 -786.5574 51.9453 -786.5574 0.0000 -786.5574
## [15] 0.0000 -786.5574 51.9453 -786.5574 0.0000 -786.5574benefit_of_controlled_burning## [1] 0.0000 6492.6617 0.0000 -786.5574 0.0000 -786.5574 0.0000
## [8] -786.5574 0.0000 -786.5574 7331.1643 -838.5026 0.0000 -786.5574
## [15] 0.0000 6492.6617 0.0000 -786.5574 0.0000 -786.5574The last step in this decision model is to compute the total benefits by summing up the yearly values. A slight complication in this is that very few people are indifferent about the time at which they are faced with costs or benefits. Most people would prefer to have income now over having that same income in the future, and they feel the opposite way about expenses (most prefer expenses later). This characteristic is known as the time preference. We also say that people ‘discount’ future costs and benefits. This is considered when calculating the so-called Net Present Value (NPV), which discounts future costs or benefits by an investor-specific discount rate. Usually, the further in the future these costs and benefits occur, the lower their perceived value. decisionSupport provides a discounting function called discount(), which implements the standard equation used by economists to compute the NPV. For example, a discount rate of 4% applied to a regular income stream of 100 dollars per year over 10 years would produce:
discount(c(rep(100,10)),discount_rate = 4)## [1] 100.00000 96.15385 92.45562 88.89964 85.48042 82.19271 79.03145
## [8] 75.99178 73.06902 70.25867If we add the additional argument calculate_NPV=TRUE, to our code, the function automatically sums up these values to produce the Net Present Value:
discount(c(rep(100,10)),discount_rate=4,calculate_NPV = TRUE)## [1] 843.5332We see that the perceived net benefit, the Net Present Value, is not 1000 dollars, as we would get by just adding up the yearly values, but only 844 dollars.
Applied to our example, and using the discount_rate variable (from the input table), this produces:
NPV_controlled_burning <- discount(benefit_of_controlled_burning,
discount_rate = discount_rate,
calculate_NPV = TRUE)With this metric, we can compare the two forest management options.
Here’s an illustration of the basic steps of the simulation: