Index.Rmd

---
title:  Decision Analysis and Forecasting for Agricultural Development
  ![](images/decisionSupport_hex_sticker.png){width=1in}  
author: "Cory Whitney, Eike Luedeling et al."
github: "CWWhitney"
logo: "images/decisionSupport_hex_sticker.png"
bibliography: 
  - bib/references.bib
  - bib/packages.bib
output: 
  learnr::tutorial:
runtime: shiny_prerendered
---

![](images/decisionSupport_hex_sticker.png){width=1in}

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
# library(devtools)
# install.packages(c("bayesplot", "DiagrammeR", "ggcorrplot", "ggplot2", "ggrepel", "ggridges", "ggthemes", "igraph", "knitr", "learnr", "pander", "plyr", "rmarkdown", "tidyverse"))
# devtools::install_github("eikeluedeling/decisionSupport")
# library(brms)
library(bayesplot)
library(decisionSupport)
library(DiagrammeR)
library(ggcorrplot)
library(ggplot2)
library(ggrepel)
library(ggridges)
library(ggthemes)
library(igraph)
library(knitr)
library(learnr)
library(pander)
library(plyr)
library(rmarkdown)
library(tidyverse)

###### make variables function #######
# make variables
make_variables <- function(est,n=1)
{ x<-random(rho=est, n=n)
for(i in colnames(x)) assign(i,
                             as.numeric(x[1,i]),envir=.GlobalEnv)
}

# make_variables(as.estimate(input_estimates))

###### Needed Data ########
# load the data (use write.csv for data created in RMD)
stakeholder <- read.csv("data/stakeholder.csv")

example_decision_inputs <- read.csv("data/example_decision_inputs.csv")

# Models ##########

###### Chile function #######

input_estimates <- read.csv("data/input_estimates.csv")


# create the model function
model_function <- function(){
  income <- Yield * Market_price
  overall_costs <- Labor_cost + Management_cost
  final_result <- income - overall_costs
  return(list(final_result = final_result))
}

# Run the Monte Carlo simulation using the model function ####
example_mc_simulation <- mcSimulation(estimate = as.estimate(input_estimates),
                                      model_function = model_function,
                                      numberOfModelRuns = 100,
                                      functionSyntax = "plainNames")

###### HAIL Risk #######

hail_estimates <- read.csv("data/hail_estimates.csv")

hail_function <- function(){
  
  yields <- vv(var_mean = yield,
               var_CV = var_CV,
               n = 20)
  
  prices <- vv(var_mean = price,
               var_CV = var_CV,
               n = 20)
  
  # use 'rep' to simulate the investment
  # only in the first year (assuming the net lasts 20 years)
  invest_costs <- c(initial_investment, rep(0, 19))
  
  # use 'chance_event' to adjust yield for potential hail
  hail_adjusted_yield <- chance_event(chance = p_hail,
                                      value_if = 0,
                                      value_if_not = yield,
                                      n = 20)
  
  # calculate profit without net
  profit_no_net <- hail_adjusted_yield*prices
  
  # calculate profit with the net
  profit_with_net <- (yields*prices)-invest_costs
  
  # use 'discount' to calculate net present value
  # 'discount_rate' is expressed in percent
  NPV_no_net <- discount(profit_no_net, discount_rate = 5, calculate_NPV = TRUE)
  NPV_net <- discount(profit_with_net, discount_rate = 5, calculate_NPV = TRUE)
  
  # calculate the overall NPV of the decision (do - don't do)
  NPV_decision <- NPV_net-NPV_no_net
  
  return(list(NPV_no_net =  NPV_no_net,
              NPV_net =  NPV_net,
              NPV_decision = NPV_decision))
}

###### Run HAIL Risk model #######
hail_mc_simulation <- mcSimulation(estimate = as.estimate(hail_estimates),
                                   model_function = hail_function,
                                   numberOfModelRuns = 100,
                                   functionSyntax = "plainNames")

######  HAIL Risk evpi #######

mcSimulation_table_hail <- data.frame(hail_mc_simulation$x,
                                      hail_mc_simulation$y[3])

evpi_hail <- multi_EVPI(mc = mcSimulation_table_hail,
                        first_out_var = "NPV_decision")


#### Apple Sheep Model ####

# Inputs

input_estimates <- data.frame(
  variable = c(
    "sheep_income",
    "sheep_cost",
    "apple_income",
    "apple_cost",
    "discount_rate"
  ),
  lower = c(3000, 1000, 30000, 15000, 10),
  median = NA,
  upper = c(5000, 2500, 60000, 30000, 10),
  distribution = c("posnorm", "posnorm", "posnorm", "posnorm", "const"),
  label = c(
    "Income from sheep (euro/year)",
    "Cost of sheep (euro/year)",
    "Income from apple (euro/year)",
    "Cost of apple (euro/year)",
    "Discount Rate"
  )
)


###### Apple Sheep #########
# Model function 

model_function <- function() {
  # Estimate the income in a normal season
  AF_income <- sheep_income + apple_income
  AF_cost <- sheep_cost + apple_cost
  # Estimate the final results from the model
  AF_final_result <- AF_income - AF_cost
  sheep_only <- sheep_income - sheep_cost
  
  Decision_benefit <- AF_final_result - sheep_only
  #Calculating NPV
  #AF System
  AF_NPV <- discount(AF_final_result,
                     discount_rate = discount_rate,
                     calculate_NPV = TRUE)#NVP of AF system
  NPV_sheep_only <- discount(sheep_only,
                             discount_rate = discount_rate,
                             calculate_NPV = TRUE) #NVP of grassland
  NPV_decision <- discount(Decision_benefit,
                           discount_rate = discount_rate,
                           calculate_NPV = TRUE)
  # Generate the list of outputs from the Monte Carlo simulation
  return(
    list(
      NPV_Agroforestry_System = AF_NPV,
      NPV_Treeless_System = NPV_sheep_only,
      NPV_decision = NPV_decision
    )
  )
}

# simulation result
apple_sheep_mc_simulation <- mcSimulation(
  estimate = as.estimate(input_estimates),
  model_function = model_function,
  numberOfModelRuns = 80,
  functionSyntax = "plainNames"
)

###### burkina function #######

input_table <- read.csv("data/example_input_table.csv")


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-calculate 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 rain-fed 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))
}
###### run burkina model #######
# second simulation
mcSimulation_results <- decisionSupport::mcSimulation(
  estimate = decisionSupport::estimate_read_csv("data/example_input_table.csv"),
  model_function = example_decision_function,
  numberOfModelRuns = 100,
  functionSyntax = "plainNames"
)

###### burkina evpi #####
mcSimulation_table <- data.frame(mcSimulation_results$x, mcSimulation_results$y[1:3])

evpi <- multi_EVPI(mc = mcSimulation_table, first_out_var = "Interv_NPV")
###### Burkina pls ########
pls_result <- plsr.mcSimulation(object = mcSimulation_results,
                                resultName = names(mcSimulation_results$y)[3], ncomp = 1)

# suppress warnings
options(warn = - 1)

# brms stuff with Schiffers
# library(brms, warn.conflicts=F, quietly=T)
# priors <- c(set_prior("normal(200, 5)", coef = "rasseN"), set_prior("normal(210, 5)", coef = "rasseK"))
# rasse <- c("K", "K", "K", "N", "N", "N", "N", "N", "N")
# ll <- c(225, 222, 221, 219, 219, 216, 221, 218, 217)
# daten <- data.frame(rasse=rasse, ll=ll)
# legeleistung <- brms::brm(formula = ll ~ rasse -1, data = daten , prior = priors, silent = 2)
```

```{r packages, warning=FALSE, include = FALSE}
#Automatically write R package citation entries to a .bib file
knitr::write_bib(c(.packages(), 
                  # 'brms',
                   'bayesplot',
                   'decisionSupport', 
                   'DiagrammeR',
                   'ggcorrplot',
                   'ggplot2',
                   'ggrepel',
                   'ggridges',
                   'ggthemes',
                   'igraph',
                   'knitr',
                   'learnr',
                   'plyr',
                   'rmarkdown',
                   'tidyverse'), 'bib/packages.bib')
```

<!-- List of Lecture RMD files to be compiled -->

```{r child="Lecture_01_Introduction.Rmd"}
```

<!-- ```{r child="Tools_Experiential.Rmd"} -->
<!-- ``` -->


```{r child="Lecture_Schedule.Rmd"}
```

```{r child="Seminar_Schedule.Rmd"}
```

```{r child="Lecture_02_Decision_Analysis.Rmd"}
```

```{r child="Lecture_03_Define_Decision.Rmd"}
```

```{r child="Lecture_04_Case_Studies.Rmd"}
```

```{r child="Lecture_05_Participatory_Methods.Rmd"}
```

```{r child="Lecture_06_Decision_Models.Rmd"}
```

```{r child="Lecture_07_Forecasts.Rmd"}
```

```{r child="Lecture_08_Practical_Schiffers.Rmd"}
```

```{r child="Lecture_08_01_Bayesian_Thinking.Rmd"}
```

```{r child="Lecture_09_calibration_lecture.Rmd"}
```


```{r child="Lecture_10_Analyst_Profile.Rmd"}
```

```{r child="Seminar_01_R_RStudio_1.Rmd"}
```

```{r child="Seminar_02_R_RStudio_2.Rmd"}
```

```{r child="Seminar_03_Model_Updates.Rmd"}
```


```{r child="Seminar_04_RMarkdown.Rmd"}
```

```{r child="Seminar_05_Git_Github.Rmd"}
```

```{r child="Seminar_06_Model_Programming.Rmd"}
```

```{r child="Seminar_07_Model_Programming_2.Rmd"}
```

```{r child="Seminar_08_Functions.Rmd"}
```

<!-- # ```{r child="Seminar_07_calibration_seminar.Rmd"} -->
<!-- # ``` -->

```{r child="Seminar_09_Value_of_information.Rmd"}
```

```{r child="Seminar_09_01_VoI.Rmd"}
```

```{r child="Seminar_10_Communicating.Rmd"}
```

<!-- ```{r child="Seminar_08_models_seminar.Rmd"} -->
<!-- ``` -->

```{r child="Seminar_11_forecast_seminar.Rmd"}
```

<!-- ```{r child="Seminar_12_Writing.Rmd"} -->
<!-- ``` -->

```{r child="Seminar_12_1_Citations.Rmd"}
```

```{r child="Impressum.Rmd"}
```

```{r child="References.Rmd"}
```