Tutorial.Rmd

---
title: "Decision Modeling for Agroecology"
author: "Cory Whitney"
date: "`r Sys.Date()`"
output: learnr::tutorial
runtime: shiny_prerendered
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(decisionSupport)
library(ggplot2)
library(learnr)
library(shiny)
library(DiagrammeR)

source("functions/monte_carlo.R")
source("functions/apples_sheep_in_class.R")
source("functions/apple_sheep_pareto.R")

```


## **Introduction**

Welcome to the interactive tutorial on **Decision Modeling for Agroecology & Conservation**. 
This tutorial will guide you through:

- **Theoretical foundations** of decision modeling 
- **Identifying decision options and outcomes** 
- **Structuring causal models** using Directed Acyclic Graphs (DAGs)
- **Monte Carlo simulation** for decision-making
- **Expert-elicited probabilities** and model validation
- **Evaluating models** with Pareto-optimality

```{r question0, echo=FALSE}
textInput("intro_question", "What do you already know about decision modeling?", placeholder="Write a short response...")
```

---

## **1. Decision Modeling**

Decision-making in agroecology often requires navigating trade-offs, for example between productivity and biodiversity. Decisions are further complicated by system complexity and uncertainty, i.e. the impacts of climate variability or market fluctuations. In this section, we will explore why decision modeling matters for agroecology, how it can guide sustainable land management, and what tools can help us make informed decisions.

```{r question1, echo=FALSE}
quiz(
  question("What is a good reason to use a probabilistic approach to support decision-making in agroecology?",
    answer("It captures uncertainty in outcomes and trade-offs.", correct = TRUE),
    answer("It provides fixed deterministic solutions.", correct = FALSE),
    answer("It replaces the need for expert knowledge.", correct = FALSE),
    answer("It avoids the need for validation.", correct = FALSE)
  )
)
```

---

## **2. Identifying Decision Options and Outcomes**

**Activity:** Identify and list **decision options** and **outcomes of interest** from the perspective of decision-makers.

```{r decision_options, echo=FALSE}
textInput("decision_options", "List the decision options relevant to an agroecological problem:", placeholder="Example: Adopting agroforestry, expanding irrigation, switching to organic inputs...")
```

```{r decision_outcomes, echo=FALSE}
textInput("decision_outcomes", "List the key outcomes of interest for decision-makers:", placeholder="Example: Yield improvement, biodiversity conservation, financial returns...")
```

---

## **3. Structuring a Causal Decision Model (DAGs)**

**Activity:** Sketch a decision model using the interactive diagram tool or upload a hand-drawn image.

### **Option 1: Interactive DAG Drawing**

Share a link of a graph you build in [Loopy](https://ncase.me/loopy/) by [Explorable Explanations](http://explorableexplanations.com/)

### **Option 2: Upload a Hand-Drawn Model**

Sketch out a decision model on a piece of paper and share with us.

```{r upload_model, echo=FALSE}
fileInput("model_upload", "Upload your decision model (optional)", accept = c("image/png", "image/jpeg"))
```

Or you can adjust this one if you are familiar with DOT language. 

```{r decision_model, exercise=TRUE}
DiagrammeR::grViz("
digraph decision_model {
  # Define node styles
  node [shape=ellipse, style=filled, fillcolor=lightyellow] Apple_Income Sheep_Income Total_incomes ;
  node [shape=ellipse, style=filled, fillcolor=lightpink] Sheep_cost Apple_cost Total_costs;
  node [shape=ellipse, style=filled, fillcolor=lightblue] Outcome_Profit;
  
  # Add layout
  graph [rankdir=LR]
  
  # Add labels
  Apple_Income [label = 'Apple Income (€)']
  Sheep_Income [label = 'Sheep Income (€)']
  Total_incomes [label = 'Total Income']
  Apple_cost [label = 'Apple Cost (€)']
  Sheep_cost [label = 'Sheep Cost (€)']
  Total_costs [label = 'Total Cost']
  Outcome_Profit [label = 'Net Profit (€)']

  # Define relationships
  Sheep_cost -> Total_costs;
  Apple_cost -> Total_costs;
  Apple_Income -> Total_incomes;
  Sheep_Income -> Total_incomes;
  Total_costs -> Outcome_Profit;
  Total_incomes -> Outcome_Profit
}
")

```

---

## **4. Implementing a Monte Carlo Model**

Now for a very simple **Monte Carlo Simulation** in R. We will analyze a **decision**. The choice of a farmer to start agroforestry (apple trees + sheep) or maintain a treeless sheep system. Use `set.seed` and change the number of simulations to 1000. 

```{r apples_model, exercise=TRUE}

# set.seed(42)
num_simulations <- 100

# Define ranges per ha per year Euro
# lower and upper
apple_income <- runif(n = num_simulations, 
                      min = 3000, 
                      max = 60000)

apple_costs <- runif(n = num_simulations, 
                      min = 15000, 
                     max = 30000)

apple_profits <- apple_income - apple_costs
```

```{r apples_model-solution}

set.seed(42)
num_simulations <- 1000

# Define ranges per ha per year Euro
# lower and upper
apple_income <- runif(n = num_simulations, 
                      min = 3000, 
                      max = 60000)

apple_costs <- runif(n = num_simulations, 
                      min = 15000, 
                     max = 30000)

apple_profits <- apple_income - apple_costs
```

Add sheep to the horticulture system. So we can judge the relative benefits of transition to a silvopastoral agroforestry system. Add the sheep and apple profits as a last step in the function.

```{r sheep_model, exercise=TRUE}

# Euro per ha per year
sheep_income <- runif(n = num_simulations, 
                      min = 2000, 
                      max = 5000)

sheep_costs <- runif(n = num_simulations, 
                     min = 1000, 
                     max = 2500)

sheep_profits <- sheep_income - sheep_costs
```

```{r sheep_model-solution}

# Euro per ha per year
sheep_income <- runif(n = num_simulations, 
                      min = 2000, 
                      max = 5000)

sheep_costs <- runif(n = num_simulations, 
                     min = 1000, 
                     max = 2500)

sheep_profits <- sheep_income - sheep_costs

total_profits <- apple_profits + sheep_profits
```

Plot a histogram of the profits. Change the color of the apple_profits to grey. 

```{r plot_profits_overlay, exercise=TRUE}
hist(total_profits, col = "white")
hist(apple_profits, add = TRUE, col = "pink")
hist(sheep_profits, add = TRUE, col = "black")
```

```{r plot_profits_overlay-solution}
hist(total_profits, col = "white")
hist(apple_profits, add = TRUE, col = "grey")
hist(sheep_profits, add = TRUE, col = "black")
```

```{r question2, echo=FALSE}
quiz(
  question("What does the histogram tell us about uncertainty?",
    answer("There is high variation in profits, suggesting high uncertainty.", correct = TRUE),
    answer("All outcomes are nearly the same, meaning low uncertainty.", correct = FALSE)
  )
)
```

---

## **5. Pareto-Optimality**

Now we will calcualte Pareto-Optimality. We will use the concept of the Pareto front to find the best trade-offs between two competing objectives. We can learn which solutions are *optimal* (better in at least one way, without being worse in another) and which solutions are *dominated* (we could improve at least one objective without making the other worse).

First we need to increase the model complexity a little by adding management options. We add stocking density (sheep per hectare), tree density (trees per hectare), proportion of tree area, and proportion of sheep area (1 - proportion_tree_area). *Adjust these numbers to reflect your best estimate.*

```{r management_variables, exercise=TRUE}

# Define management variables per hectare
max_stocking_density <- 20  # Max sheep per ha before overgrazing effects occur
stocking_density <- runif(num_simulations, min = 5, max = max_stocking_density)

max_tree_density <- 500   # Maximum trees per hectare
tree_density <- runif(num_simulations, min = 50, max = max_tree_density)

# Define proportion of tree area vs. sheep area
proportion_tree_area <- runif(num_simulations, min = 0.2, max = 0.8)
proportion_sheep_area <- 1 - proportion_tree_area  # Remaining area for sheep
```

```{r management_variables-solution}

# Define management variables per hectare
max_stocking_density <- 40  # Max sheep per ha before overgrazing effects occur
stocking_density <- runif(num_simulations, min = 5, max = max_stocking_density)

max_tree_density <- 200   # Maximum trees per hectare
tree_density <- runif(num_simulations, min = 50, max = max_tree_density)

# Define proportion of tree area vs. sheep area
proportion_tree_area <- runif(num_simulations, min = 0.2, max = 0.8)
proportion_sheep_area <- 1 - proportion_tree_area  # Remaining area for sheep

total_profits_agroforestry <- apple_profits + (sheep_profits * proportion_sheep_area)
total_profits_treeless <- sheep_profits  # No tree component in treeless system
```

Adjust these biodiversity numbers to reflect your best estimate. For example, depending on the ecosystem, each tree likely adds far fewer species and each sheep likely has far less of a negative impact. 

```{r biodiversity_variables, exercise=TRUE}

# Define species richness and grazing impact factors per hectare
species_richness_per_tree <- 0.05  # Each tree adds 0.05 species per ha
species_loss_per_sheep <- 0.02  # Each sheep reduces 0.02 species per ha

# Calculate biodiversity per ha
tree_species_richness <- tree_density * species_richness_per_tree
grazing_species_loss <- stocking_density * species_loss_per_sheep

# Compute biodiversity outcomes per ha
biodiversity_agroforestry <- tree_species_richness - grazing_species_loss
biodiversity_treeless <- -grazing_species_loss  # No trees, only grazing effect

# Calculate difference (decision impact of trees)
# measures the financial difference when switching to trees
total_profits <- total_profits_agroforestry - total_profits_treeless

# Calculate difference (decision impact of trees on biodiversity)
# measures the net gain/loss in biodiversity from switching to trees
biodiversity_impact <- biodiversity_agroforestry - biodiversity_treeless

# Define a weighting factor for biodiversity impact (convert biodiversity into € equivalent)
biodiversity_value_per_species <- 500  # Assume each additional species is worth €500

```

```{r biodiversity_variables-solution}

# Define species richness and grazing impact factors per hectare
species_richness_per_tree <- 0.01  # Each tree adds 0.001 species per ha
species_loss_per_sheep <- 0.001  # Each sheep reduces 0.001 species per ha

# Calculate biodiversity per ha
tree_species_richness <- tree_density * species_richness_per_tree
grazing_species_loss <- stocking_density * species_loss_per_sheep

# Compute biodiversity outcomes per ha
biodiversity_agroforestry <- tree_species_richness - grazing_species_loss
biodiversity_treeless <- -grazing_species_loss  # No trees, only grazing effect

# Calculate difference (decision impact of trees)
# measures the financial difference when switching to trees
total_profits <- total_profits_agroforestry - total_profits_treeless

# Calculate difference (decision impact of trees on biodiversity)
# measures the net gain/loss in biodiversity from switching to trees
biodiversity_impact <- biodiversity_agroforestry - biodiversity_treeless

# Define a weighting factor for biodiversity impact (convert biodiversity into € equivalent)
biodiversity_value_per_species <- 500  # Assume each additional species is worth €500


```

For the Pareto front plot we will first add the data together, store results in dataframe and plot a scatter plot. Change the color of the points in the scatter plot to green. 

```{r results_data, exercise=TRUE}

results <- data.frame(
  total_profits,
  biodiversity_impact
)

# Scatter plot
plot(results$biodiversity_impact, results$total_profits,
     col = "blue", pch = 19, xlab = "Biodiversity Impact (Species Change per ha)",
     ylab = "Total Profits (€ per ha)", main = "Pareto Front: Agroforestry vs. Treeless System")
```

```{r results_data-solution}

results <- data.frame(
  total_profits,
  biodiversity_impact
)

# Scatter plot
plot(results$biodiversity_impact, results$total_profits,
     col = "green", pch = 19, xlab = "Biodiversity Impact (Species Change per ha)",
     ylab = "Total Profits (€ per ha)", main = "Pareto Front: Agroforestry vs. Treeless System")
```

Now that we see what the scatter plot looks like what is your impression? 

```{r question3, echo=FALSE}
textInput("pareto_question", "What do you expect the trade-offs to be between profit and biodiversity?", 
          placeholder="Write your prediction before running the analysis...")
```

Now we can look for the non-dominated points where no other solution is better in both biodiversity and total profits. A point (biodiversity, profit) is Pareto-optimal if no other point has higher biodiversity and higher profits. Change the color of the points to blue.

```{r pareto_points, exercise=TRUE}

# Sort results by biodiversity impact
results <- results[order(results$biodiversity_impact, decreasing = TRUE),]

# Initialize Pareto front
pareto_front <- results[1, ]  # Start with the first point

# Loop through sorted results and keep only non-dominated points
for (i in 2:nrow(results)) {
  if (results$total_profits[i] > tail(pareto_front$total_profits, 1)) {
    pareto_front <- rbind(pareto_front, results[i, ])
  }
}

# Plot all points
plot(results$biodiversity_impact, results$total_profits,
     col = "green", pch = 19, xlab = "Biodiversity Impact (Species Change per ha)",
     ylab = "Total Profits (€ per ha)", main = "Pareto Front: Agroforestry vs. Treeless System")

# Highlight Pareto front points
points(pareto_front$biodiversity_impact, pareto_front$total_profits,
       col = "red", pch = 19, type = "o", lwd = 2)
```

```{r pareto_points-solution}

# Sort results by biodiversity impact
results <- results[order(results$biodiversity_impact, decreasing = TRUE),]

# Initialize Pareto front
pareto_front <- results[1, ]  # Start with the first point

# Loop through sorted results and keep only non-dominated points
for (i in 2:nrow(results)) {
  if (results$total_profits[i] > tail(pareto_front$total_profits, 1)) {
    pareto_front <- rbind(pareto_front, results[i, ])
  }
}

# Plot all points
plot(results$biodiversity_impact, results$total_profits,
     col = "green", pch = 19, xlab = "Biodiversity Impact (Species Change per ha)",
     ylab = "Total Profits (€ per ha)", main = "Pareto Front: Agroforestry vs. Treeless System")

# Highlight Pareto front points
points(pareto_front$biodiversity_impact, pareto_front$total_profits,
       col = "blue", pch = 19, type = "o", lwd = 2)
```

Now we can zoom in on the Pareto-optimal points and examine the conditions for key operational variables we calculated earlier.

```{r pareto_conditions, exercise=TRUE}

# Extract operational conditions for Pareto-optimal points
pareto_conditions <- results[results$biodiversity_impact %in% pareto_front$biodiversity_impact &
                               results$total_profits %in% pareto_front$total_profits,]

# Merge with original dataset to include operational variables
pareto_conditions <- merge(pareto_conditions, data.frame(
  biodiversity_impact, total_profits, stocking_density, tree_density, proportion_tree_area),
  by = c("biodiversity_impact", "total_profits"))

# Show Pareto-optimal conditions
print(pareto_conditions)
```

```{r pareto_conditions-solution}

# Extract operational conditions for Pareto-optimal points
pareto_conditions <- results[results$biodiversity_impact %in% pareto_front$biodiversity_impact &
                               results$total_profits %in% pareto_front$total_profits,]

# Merge with original dataset to include operational variables
pareto_conditions <- merge(pareto_conditions, data.frame(
  biodiversity_impact, total_profits, stocking_density, tree_density, proportion_tree_area),
  by = c("biodiversity_impact", "total_profits"))

# Show Pareto-optimal conditions
print(pareto_conditions)
```

Finally we can plot the Pareto. Make the color of the points blue.

```{r plot_pareto, exercise=TRUE}
# Determine dynamic position based on plot coordinates
text_positions <- ifelse(
  pareto_front$biodiversity_impact > median(pareto_front$biodiversity_impact), 2, 4
)

# Adjust for very high or low points
text_positions <- ifelse(
  pareto_front$total_profits > quantile(pareto_front$total_profits, 0.75), 1, text_positions
)
text_positions <- ifelse(
  pareto_front$total_profits < quantile(pareto_front$total_profits, 0.25), 3, text_positions
)

# Plot Pareto front
plot(pareto_front$biodiversity_impact, pareto_front$total_profits,
     col = "red", pch = 19, type = "o", lwd = 2,
     xlab = "Biodiversity Impact (Species Change per ha)",
     ylab = "Total Profits (€ per ha)", main = "Pareto Front: Agroforestry vs. Treeless System")

# Add footnote at the bottom
mtext("S = Stocking Density (sheep/ha), T = Tree Density (trees/ha), P = Proportion of Tree Area", 
      side = 1, line = 4, cex = 0.8)

# Create multi-line labels with line breaks
labels <- paste0("S:", round(pareto_conditions$stocking_density, 1), "\n",
                 "T:", round(pareto_conditions$tree_density, 1), "\n",
                 "P:", round(pareto_conditions$proportion_tree_area, 2))

# Annotate with dynamically positioned text
text(pareto_front$biodiversity_impact, pareto_front$total_profits,
     labels = labels, pos = text_positions, cex = 0.8)
```

```{r plot_pareto-solution}
# Determine dynamic position based on plot coordinates
text_positions <- ifelse(
  pareto_front$biodiversity_impact > median(pareto_front$biodiversity_impact), 2, 4
)

# Adjust for very high or low points
text_positions <- ifelse(
  pareto_front$total_profits > quantile(pareto_front$total_profits, 0.75), 1, text_positions
)
text_positions <- ifelse(
  pareto_front$total_profits < quantile(pareto_front$total_profits, 0.25), 3, text_positions
)

# Plot Pareto front
plot(pareto_front$biodiversity_impact, pareto_front$total_profits,
     col = "blue", pch = 19, type = "o", lwd = 2,
     xlab = "Biodiversity Impact (Species Change per ha)",
     ylab = "Total Profits (€ per ha)", main = "Pareto Front: Agroforestry vs. Treeless System")

# Add footnote at the bottom
mtext("S = Stocking Density (sheep/ha), T = Tree Density (trees/ha), P = Proportion of Tree Area", 
      side = 1, line = 4, cex = 0.8)

# Create multi-line labels with line breaks
labels <- paste0("S:", round(pareto_conditions$stocking_density, 1), "\n",
                 "T:", round(pareto_conditions$tree_density, 1), "\n",
                 "P:", round(pareto_conditions$proportion_tree_area, 2))

# Annotate with dynamically positioned text
text(pareto_front$biodiversity_impact, pareto_front$total_profits,
     labels = labels, pos = text_positions, cex = 0.8)
```

```{r question4, echo=FALSE}
textInput("pareto_reflection", "What did you observe in the Pareto front?")
```

---

## **6. Next Steps**

```{r own_work, exercise=TRUE}
# Define a decision problem related to agroecology. 
# What are the key variables? 

# Example: Define key income and cost variables
your_income <- runif(100, min = ???, max = ???)  # Fill in realistic values
your_costs <- runif(100, min = ???, max = ???) 

# Compute profits
your_profits <- your_income - your_costs

# --- Next Steps ---
# Add an uncertainty factor (e.g., climate risk)
# Include biodiversity or sustainability impacts
# run a Monte Carlo simulation
```

---

## **7. Final Reflection**

- What insights did you gain from this tutorial?
- How would you improve your decision model?
- How does this relate to **real-world agroecology and conservation?**

```{r feedback1, echo=FALSE}
sliderInput("engagement", "How engaging was this tutorial?", min = 1, max = 5, value = 3)
```

```{r feedback2, echo=FALSE}
textInput("difficulty", "What was the hardest part to understand?")
```

```{r feedback3, echo=FALSE}
textInput("feedback", "Share your key takeaways:", placeholder="Write your reflections here...")
```

---

## **Next Steps**

- Apply Decision Analysis to your own research
- Continue refining expert-elicited models

Thanks for following this tutorial. As thanks, here is an open space for running your scripts. 

```{r blank, exercise=TRUE}

```

---

## References

Whitney, C, K Shepherd, and E Luedeling. “Decision Analysis Methods Guide; Agricultural Policy for Nutrition.” World Agroforestry (ICRAF) Working Paper series, no. 275 (2018): 40. http://dx.doi.org/10.5716/WP18001.PDF.

Whitney, Cory, Denis Lanzanova, Caroline Muchiri, Keith D. Shepherd, Todd S. Rosenstock, Michael Krawinkel, John R. S. Tabuti, and Eike Luedeling. “Probabilistic Decision Tools for Determining Impacts of Agricultural Development Policy on Household Nutrition.” Earth’s Future 6, no. 3 (2018): 359–72. https://doi.org/10.1002/2017EF000765.

Whitney, Cory, Gordon O’Brien, Vuyisile Dlamini, Ikhothatseng Jacob Greffiths, Chris Dickens, and Eike Luedeling. “Balancing Ecosystem Sustainability and Irrigated Smallholder Agriculture: A Modeling Approach for Water Resource Management.” Journal of Hydrology 651 (April 2025): 132560. https://doi.org/10.1016/j.jhydrol.2024.132560.

Whitney, Cory, Lisa Biber-Freudenberger, and Eike Luedeling. “Decision Analytical Methods for Assessing the Efficacy of Agroecology Interventions.” CABI Agriculture and Bioscience 4, no. 1 (2023): 9. https://doi.org/10.1186/s43170-023-00151-9.