This document illustrates the use of various functions within the decisionSupport R package (Luedeling, Goehring, and Schiffers 2019). We use these to develop a cost-benefit model and implement a Monte Carlo simulation with the example of agroforestry interventions promoted in Northwest Vietnam (see Do et al. 2020). We generate a model to calculate the Net Present Value of proposed agroforestry systems as compared to the common conventional monoculture cropping systems of maize.

We use the decisonSupport() function to perform probabilistic Monte Carlo simulations. The function performs calculations using distributions and ranges which is a major advantage over traditional spreadsheet calculations. In a Monte Carlo simulation, a particular calculation is repeated many times. Each time, a random draw from defined distribution of each variable are used as inputs for calculating the function. The output then consists of a set of values instead of single estimates.

The function requires two arguments:

  1. a data table (.csv format) specifying the names and distributions for all uncertain variables and,

  2. an R function that predicts decision outcomes based on the variables named in the data table.

The chance_event() function calculates the time series occurrence of particular events (risk events in this case). This function requires two main arguments:

  1. the probability that a particular risk event will happen in a given year and

  2. the time line for the simulation.

The gompertz_yield() function calculates the yields of tree crops using identified parameters given an assumed pattern.

The value varier function vv() produces realistic time series of data that includes year-to-year variation around a specified mean. This is determined by a desired coefficient of variation. Both mean and coefficient of variation are provided to the model via the input table.

The discount() function calculated NPV in the last step of the decision modeling procedure. The function implements discounting of future costs or benefits by a decision maker-specific discount rate, which is also specified in the input tables.

The cost-benefit model AF_benefit calculates the NPV:

AF_benefit <- function(x, varnames)
{
#Define each variable as vectors of 25 values corresponding to 25 years of simulation 
  mono_costs <- rep(0, n_years) #the annual cost of maize monoculture 
  establisment_cost <- rep(0, n_years) #costs for setting up agroforestry systems (AF)
  maintenance_cost <- rep(0, n_years) #the annual cost spent on AF
  tree_yield <- rep(0, n_years) #the yield of longan
  AF_maize <- rep(0, n_years) #the yield of maize in AF
  mono_maize <- rep(0, n_years) #the yield of maize in monocultures
  grass_benefit <- rep(0, n_years) #the benefit from grass
  reduced_soil_loss <- rep(0, n_years) #amount of soil loss reduced in AF compared to monoculture systems
  
 #Simulate the chance for risk events to occur during the simulation period
  tree.risk <- chance_event(tree_risk, value_if = 1, n = n_years)
  maize.risk <- chance_event(maize_risk, value_if = 1, n = n_years)
  
  #Calculate system costs and convert currency from VND to USD using currency 
  #exchange rate (cur_change) applied for any variable with monetary value 
  #(i.e. costs and benefits)
  maintenance_cost[1] <- 0
  maintenance_cost[2:4] <- vv(main_cost_stage1, CV_cost, 3)/cur_change 
  maintenance_cost[5:10] <- vv(main_cost_stage2, CV_cost, 6)/cur_change
  maintenance_cost[11:25] <- vv(main_cost_stage3, CV_cost, 15)/cur_change
  
  establisment_cost[1] <- establis_cost/cur_change
  establisment_cost[2:25] <- 0


  ##Calculate system benefits
  
  #benefit of longan fruit
  tree_yield <- gompertz_yield(max_harvest = max_tree_harvest, 
                             time_to_first_yield_estimate = time_first_tree_est, 
                             time_to_second_yield_estimate = time_sec_tree_est, 
                             first_yield_estimate_percent = first_tree_est_per, 
                             second_yield_estimate_percent = sec_tree_est_per, 
                             n_years = n_years, 
                             var_CV = CV_tree_yield, 
                             no_yield_before_first_estimate = FALSE)
  
  tree_yield <- tree_yield*(1-tree.risk*yield_tree_risk)
  tot_tree_yield <- tree_yield*num_of_tree
  tree_benefit <- (tot_tree_yield*tree_price)/cur_change
  
  #benefit of maize in AF systems
  time <- 1:n_years
  decay_speed_AF <- -log(1-decay_rate_AF)
  AF_maize <- min_AFmaize+(max_AFmaize-min_AFmaize)*exp(-decay_speed_AF*(time-1))
  tot_AF_maize <- vv(AF_maize*(1-maize.risk*yield_maize_risk), CV_maize_yield, n_years)
  AF_maize_benefit <- (tot_AF_maize*vv(maize_price, CV_maize_price, n_years))/cur_change

  
  #benefit of maize in monoculture
  decay_speed_mono <- -log(1-decay_rate_mono)
  mono_maize <- max_monomaize*exp(-decay_speed_mono*(time-1))
  tot_mono_maize <- vv(mono_maize*(1-maize.risk*yield_maize_risk), CV_maize_yield, n_years)
  mono_revenue <- (tot_mono_maize*vv(maize_price, CV_maize_price, n_years))/cur_change
  mono_costs <- vv(mono_cost, CV_cost, n_years)/cur_change
  mono_benefit <- mono_revenue-mono_costs

  
  #benefit of grass
  grass_benefit <- vv(grass_profit, CV_grass_profit, 25)/cur_change
  
  #benefit of soil erosion control
  reduced_soil_loss[1] <- 0
  reduced_soil_loss[2:10] <- reduced_percents1*mono_soil_loss
  reduced_soil_loss[10:25] <- reduced_percents2*mono_soil_loss
  soil_saved <- vv(reduced_soil_loss, CV_eros_control, n = n_years)
  erosion_benefit <- (soil_saved*soil_saved_payment)/cur_change

# Calculate NPVs of AF and monoculture and NPV of AF vs. monocultures of maize
  AF_revenue <- tree_benefit + AF_maize_benefit + grass_benefit+erosion_benefit 
  AF_benefit <- AF_revenue-maintenance_cost-establisment_cost

  cash_flow <- discount(AF_benefit, 
                      discount_rate = discount_rate, 
                      calculate_NPV = FALSE) 
  
#cummulative cash flow of agroforestry system  
  cum_cash_flow <- cumsum(cash_flow)

#NPV of monoculture system
  NPV_mono <- discount(mono_benefit, discount_rate = discount_rate, 
                       calculate_NPV = TRUE ) 

#Benefit of choosing agroforestry over maize monoculture
  tradeoff_benefit <- AF_benefit-mono_benefit

#NPV of AF system
  NPV_AF <- discount(AF_benefit, 
                   discount_rate = discount_rate, 
                   calculate_NPV = TRUE)

#final NPV of the decision to choose AF or monocultures
  NPV_tradeoff <- discount(tradeoff_benefit, 
                         discount_rate = discount_rate, 
                         calculate_NPV = TRUE )
  
#in the return list, one could indicate any outcome that they wish to see from the model
  return(list(trade_off = NPV_tradeoff, 
            NPV_maize = NPV_mono, 
            NPV = NPV_AF))
}

The probabilistic simulation was implemented using the decisionSupport function to return a Monte Carlo data table with 10, 000 values for each of the input and response variables in the return list. A VIP score generated for each input variable of the PLS regression indicate its influenctial importance contriduting to model outcome.

decisionSupport(inputFilePath = paste("Longan+maize+grass.csv", sep = ""), 
                outputPath = paste("MCResults", sep = ""), 
                write_table = TRUE, 
                welfareFunction = AF_benefit, 
                numberOfModelRuns = 10000, 
                functionSyntax = "plainNames")

The outcome table is then used to calculate the Value of Information for uncertain variables using the multi_EVPI function in the decisionSupport package (Luedeling, Goehring, and Schiffers 2019):

The multi_EVPI() function calculates the expected value of perfect information (EVPI) for multiple variables in the model using the simulation outcome table from the decisionSupport()function.

MCall <- read.table(file = "MCResults/mcSimulationResults.csv", 
                    header = TRUE, sep  = ",")

multi_EVPI(MCall, "trade_off", write_table = TRUE)

The full NPV distributions of agroforestry and maize monoculture acquired from MC simulation look like this:

VIP scores from the sensitivity analysis (left) and Value of Information analysis (right):

References

Luedeling, Eike, Lutz Goehring, and Katja Schiffers. 2019. DecisionSupport: Quantitative Support of Decision Making Under Uncertainty. http://www.worldagroforestry.org/.