Seminar_09_01_VoI.Rmd

## Seminar 9.1: Value of Information {#voi_1} 
<!-- reference with [Models](#value_of_information_seminar_part2) -->

Welcome to the second part of the Value of Information seminar for the course **Decision Analysis and Forecasting for Agricultural Development**.Feel free to bring up any questions or concerns in the Slack or to [Dr. Cory Whitney](mailto:cory.whitney@uni-bonn.de?subject=[Seminar_4]%20Decision%20Analysis%20Lecture) or the course tutor.

### Practical Value of Information example:

Let's use a modified version of sheep apple agroforestry from Seminar 7 (#model_programming_cont). Just for this class let's presume that the farmer has a grassland where they graze sheep. Intervention: implement apple agroforestry. Let's prepare our input estimates for the system.

```{r libraries_voi_part2, eval=FALSE, echo=TRUE}
library(ggplot2)
library(decisionSupport)
```


```{r input_table_voi_test, exercise=TRUE}
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"
  )
)

# show the result
input_estimates
```

```{r input_table_voi_test-solution}
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"
  )
)

# show the result
input_estimates
```

If you want you can write the input estimates into a local file on your machine. 

```{r, eval=FALSE, echo=TRUE}
# write the results to a .csv file
df <- data.frame(input_estimates)
write.csv(df, file = "input_table.csv", row.names = FALSE)
```

In the next step, let's build a function to model the system and calculate the NPV of the AF system vs only sheep. Simulate it using the `mcSimulation` function from the `decisionSupport` package with 5000 'numberOfModelRuns'. Using 'plot_distributions' plot NPV of intervention trade off vs no intervention.

```{r apple_sheep_voi_test, exercise=TRUE}
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
  # baseline with sheep only
  sheep_only <- sheep_income - sheep_cost
  # should I plant trees in the sheep pastures? 
  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
    )
  )
}
```

```{r apple_sheep_voi_test-solution}
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
    )
  )
}
apple_sheep_mc_simulation <- mcSimulation(
  estimate = as.estimate(input_estimates),
  model_function = model_function,
  numberOfModelRuns = 80,
  functionSyntax = "plainNames"
)

plot_distributions(
  mcSimulation_object = apple_sheep_mc_simulation,
  vars = c("NPV_decision", "NPV_Treeless_System"),
  method = 'smooth_simple_overlay',
  base_size = 7,
  x_axis_name = "Outcome of AF intervention",
  scale_x_continuous(
    labels = function(x)
      x / 100000
  ),
  ggtitle("Net Present Value of the system"),
  legend.position = "bottom"
)
```

Now, let's perform sensitivity analysis using Partial Least Squares (PLS) regression using 'plsr' function. Plot the results showing a `cut_off_line` at 1 and `threshold` of 0.5. 

```{r apple_sheep_pls_test, exercise=TRUE}

pls_result_AF <- plsr.mcSimulation(
  object = apple_sheep_mc_simulation,
  resultName = names(apple_sheep_mc_simulation$y)[3],
  ncomp = 1
)

```

```{r apple_sheep_pls_test-solution}

pls_result_AF <- plsr.mcSimulation(
  object = apple_sheep_mc_simulation,
  resultName = names(apple_sheep_mc_simulation$y)[3],
  ncomp = 1
)

plot_pls(pls_result_AF,
  input_table = input_estimates,
  cut_off_line = 1,
  threshold = 0.5)

```

Save it on your machine with `ggsave`.

```{r, eval=FALSE, echo=TRUE}

plot_pls(pls_result_AF,
  input_table = input_estimates,
  cut_off_line = 1,
  threshold = 0.5)

ggsave(filename = "PLS_sheep_vs_AppleAF.png",
  plot = last_plot(),
  width = 5, 
  height = 3)

```

Now, let's quantify Value of Information by calculating expected value of perfect information (EVPI). Save the results of your MC simulation to a data frame. 

```{r assign_results_voi, exercise=TRUE}
# assign the results of x and y to a data frame

df <- data.frame(apple_sheep_mc_simulation$x, apple_sheep_mc_simulation$y[1:3])

```

```{r assign_results_voi-solution}
# assign the results of x and y to a data frame

df <- data.frame(apple_sheep_mc_simulation$x, apple_sheep_mc_simulation$y[1:3])

```

Generate the Variable Importance in the Projection (VIP) scores of the PLS for the input variables based on the apple sheep model.

```{r apple_sheep_evpi_test, exercise=TRUE}

# run evpi on the NPVs (the last few in the table, starting with "NPV_decision")

EVPI <- multi_EVPI(mc = df, first_out_var = "NPV_decision")

# plot the EVPI results for the decision
plot_evpi(EVPIresults = EVPI, decision_vars = "NPV_decision")

```

```{r apple_sheep_evpi_test-solution}

# run evpi on the NPVs (the last few in the table, starting with "NPV_decision")

EVPI <- multi_EVPI(mc = df, first_out_var = "NPV_decision")

# plot the EVPI results for the decision
plot_evpi(EVPIresults = EVPI, decision_vars = "NPV_decision")

```

Share your result with us in Slack and explain the result.