## ----include = F-------------------------------------------------------------- #set global options for knitr chunks knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width=5, fig.height=3.5 ) ## ----warning = F, include = F------------------------------------------------- #Automatically write R package citation entries to a .bib file knitr::write_bib(c(.packages(), 'chillR', 'dplyr', 'patchwork', 'plyr', 'tidyverse', 'ggplot2', 'decisionSupport'), 'packages.bib') ## ----eval = FALSE------------------------------------------------------------- # install.packages("decisionSupport") ## ----------------------------------------------------------------------------- library(decisionSupport) ## ----model-------------------------------------------------------------------- example_decision_function <- function(x, varnames){ # calculate ex-ante risks: impact the implementation of interventions #### intervention_NonPopInvolvEvent <- chance_event(intervention_NonPopInvolv, 1, 0, n = 1) # pre-calculation of common random draws for all intervention model runs #### # profits from Tropical Livestock Units (TLU) TLU <- vv(TLU_no_intervention, var_CV, n_years) TLU_profit <- vv(profit_per_TLU, var_CV, n_years) # benefits of fruit precalc_intervention_fruit_benefits <- vv(intervention_fruit_area_ha, var_CV, n_years) * vv(intervention_fruit_yield_t_ha, var_CV, n_years) * vv(intervention_fruit_profit_USD_t, var_CV, n_years) # benefits of vegetables precalc_intervention_vegetable_benefits <- vv(intervention_vegetable_area_ha, var_CV, n_years) * vv(intervention_vegetable_yield_t_ha, var_CV, n_years) * vv(intervention_vegetable_profit_USD_t, var_CV, n_years) # benefits of rainfed crops precalc_intervention_rainfed_crop_benefits <- vv(intervention_rainfed_crop_area_ha, var_CV, n_years) * vv(intervention_rainfed_crop_yield_t_ha, var_CV, n_years) * vv(intervention_rainfed_crop_profit_USD_t, var_CV, n_years) # Intervention #### for (decision_intervention_strips in c(FALSE,TRUE)) { if (decision_intervention_strips) { intervention_strips <- TRUE intervention_strips_PlanningCost <- TRUE intervention_strips_cost <- TRUE } else { intervention_strips <- FALSE intervention_strips_PlanningCost <- FALSE intervention_strips_cost <- FALSE } if (intervention_NonPopInvolvEvent) { intervention_strips <- FALSE intervention_strips_cost <- FALSE } # Costs #### if (intervention_strips_cost) { cost_intervention_strips <- intervention_adaptation_cost + intervention_tech_devices_cost + intervention_nursery_cost + intervention_wells_cost + intervention_training_cost + intervention_mngmt_oprt_cost + intervention_mngmt_follow_cost + intervention_mngmt_audit_cost } else cost_intervention_strips <- 0 if (intervention_strips_PlanningCost) { plan_cost_intervention_strips <- intervention_communication_cost + intervention_zoning_cost } else plan_cost_intervention_strips <- 0 maintenance_cost <- rep(0, n_years) if (intervention_strips) maintenance_cost <- maintenance_cost + vv(maintenance_intervention_strips, var_CV, n_years) intervention_cost <- maintenance_cost intervention_cost[1] <- intervention_cost[1] + cost_intervention_strips + plan_cost_intervention_strips # Benefits from cultivation in the intervention strips #### intervention_fruit_benefits <- as.numeric(intervention_strips) * precalc_intervention_fruit_benefits intervention_vegetable_benefits <- as.numeric(intervention_strips) * precalc_intervention_vegetable_benefits intervention_rainfed_crop_benefits <- as.numeric(intervention_strips) * precalc_intervention_rainfed_crop_benefits # Total benefits from crop production (agricultural development and riparian zone) #### crop_production <- intervention_fruit_benefits + intervention_vegetable_benefits + intervention_rainfed_crop_benefits # Benefits from livestock #### # The following allows considering that intervention strips may # restrict access to the reservoir for livestock. if (intervention_strips) TLU_intervention <- TLU * (1 + change_TLU_intervention_perc / 100) else TLU_intervention <- TLU if (decision_intervention_strips){ livestock_benefits <- TLU_intervention * TLU_profit total_benefits <- crop_production + livestock_benefits net_benefits <- total_benefits - intervention_cost result_interv <- net_benefits} if (!decision_intervention_strips){ livestock_benefits <- TLU_no_intervention * TLU_profit total_benefits <- livestock_benefits net_benefits <- total_benefits - intervention_cost result_n_interv <- net_benefits} } #close intervention loop bracket NPV_interv <- discount(result_interv, discount_rate, calculate_NPV = TRUE) NPV_n_interv <- discount(result_n_interv, discount_rate, calculate_NPV = TRUE) # Beware, if you do not name your outputs (left-hand side of the equal sign) in the return section, # the variables will be called output_1, _2, etc. return(list(Interv_NPV = NPV_interv, NO_Interv_NPV = NPV_n_interv, NPV_decision_do = NPV_interv - NPV_n_interv, Cashflow_decision_do = result_interv - result_n_interv)) } ## ----mcSimulation------------------------------------------------------------- mcSimulation_results <- decisionSupport::mcSimulation( estimate = decisionSupport::estimate_read_csv("example_input_table.csv"), model_function = example_decision_function, numberOfModelRuns = 1e3, #run 1,000 times functionSyntax = "plainNames" ) ## ----------------------------------------------------------------------------- decisionSupport::plot_distributions(mcSimulation_object = mcSimulation_results, vars = c("Interv_NPV", "NO_Interv_NPV"), method = 'smooth_simple_overlay', base_size = 7) ## ----plot_distributions_boxplot----------------------------------------------- decisionSupport::plot_distributions(mcSimulation_object = mcSimulation_results, vars = c("Interv_NPV", "NO_Interv_NPV"), method = 'boxplot') ## ----plot_distributions_box_dens---------------------------------------------- decisionSupport::plot_distributions(mcSimulation_object = mcSimulation_results, vars = "NPV_decision_do", method = 'boxplot_density') ## ----plot_cashflow------------------------------------------------------------ plot_cashflow(mcSimulation_object = mcSimulation_results, cashflow_var_name = "Cashflow_decision_do") ## ----------------------------------------------------------------------------- pls_result <- plsr.mcSimulation(object = mcSimulation_results, resultName = names(mcSimulation_results$y)[3], ncomp = 1) ## ----------------------------------------------------------------------------- input_table <- read.csv("example_input_table.csv") plot_pls(pls_result, input_table = input_table, threshold = 0) ## ----evpi, message = FALSE---------------------------------------------------- #here we subset the outputs from the mcSimulation function (y) by selecting the correct variables # this should be done by the user (be sure to run the multi_EVPI only on the variables that the user wants) mcSimulation_table <- data.frame(mcSimulation_results$x, mcSimulation_results$y[1:3]) evpi <- multi_EVPI(mc = mcSimulation_table, first_out_var = "Interv_NPV") ## ----evpi_plot---------------------------------------------------------------- plot_evpi(evpi, decision_vars = "NPV_decision_do") ## ----compound_figure, fig.width=7, fig.height=4------------------------------- compound_figure(mcSimulation_object = mcSimulation_results, input_table = input_table, plsrResults = pls_result, EVPIresults = evpi, decision_var_name = "NPV_decision_do", cashflow_var_name = "Cashflow_decision_do", base_size = 7)