---
title: "Sensitivity Analysis with Custom Efficiency Factors"
author: "Olatunde D. Akanbi"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Sensitivity Analysis with Custom Efficiency Factors}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)
```

## Overview

This vignette demonstrates how to conduct sensitivity analysis using custom nutrient efficiency factors in the manureshed package. While standard analyses use established efficiency factors (0.5 for nitrogen, 1.0 for phosphorus), users may want to explore how different efficiency assumptions affect classification results.

## Standard vs. Custom Efficiency Factors

**Standard Factors:**
- **Nitrogen: 0.5 (50%)** - Reflects typical losses during nutrient cycling, application, and uptake
- **Phosphorus: 1.0 (100%)** - Assumes full phosphorus availability

**When to Use Custom Factors:**
- Evaluating improved management practices (higher efficiency)
- Accounting for regional differences in application methods
- Conducting uncertainty analysis
- Testing different policy scenarios

## Basic Sensitivity Analysis

### Load Required Data
```{r load_data, eval=FALSE}
library(manureshed)

# Load built-in data for HUC8 watersheds
nugis_data <- load_builtin_nugis("huc8", 2016)
boundaries <- load_builtin_boundaries("huc8")
```

### Compare Different Nitrogen Efficiencies
```{r nitrogen_sensitivity, eval=FALSE}
# Standard analysis (50% efficiency)
results_standard <- agri_classify_complete_custom(
  nugis_data, 
  scale = "huc8",
  n_efficiency = 0.5,
  p_efficiency = 1.0
)

# Improved management scenario (70% efficiency)
results_high <- agri_classify_complete_custom(
  nugis_data,
  scale = "huc8", 
  n_efficiency = 0.7,
  p_efficiency = 1.0
)

# Conservative scenario (30% efficiency)
results_low <- agri_classify_complete_custom(
  nugis_data,
  scale = "huc8",
  n_efficiency = 0.3,
  p_efficiency = 1.0
)

# Compare nitrogen classifications
message("Standard (50% efficiency):\n")
print(table(results_standard$N_class))

message("\nImproved Management (70% efficiency):\n")
print(table(results_high$N_class))

message("\nConservative (30% efficiency):\n")
print(table(results_low$N_class))
```

### Visualize Sensitivity Results
```{r visualize_sensitivity, eval=FALSE}
library(ggplot2)
library(dplyr)

# Create comparison data frame
efficiency_scenarios <- data.frame(
  Efficiency = c(0.3, 0.5, 0.7),
  Scenario = c("Conservative", "Standard", "Improved"),
  Source = c(
    sum(results_low$N_class == "Source", na.rm = TRUE),
    sum(results_standard$N_class == "Source", na.rm = TRUE),
    sum(results_high$N_class == "Source", na.rm = TRUE)
  ),
  Sink_Deficit = c(
    sum(results_low$N_class == "Sink_Deficit", na.rm = TRUE),
    sum(results_standard$N_class == "Sink_Deficit", na.rm = TRUE),
    sum(results_high$N_class == "Sink_Deficit", na.rm = TRUE)
  )
)

# Plot source watersheds across scenarios
ggplot(efficiency_scenarios, aes(x = Efficiency, y = Source)) +
  geom_line(size = 1.2, color = "#2166ac") +
  geom_point(size = 3, color = "#2166ac") +
  labs(
    title = "Impact of Nitrogen Efficiency on Source Watersheds",
    x = "Nitrogen Efficiency Factor",
    y = "Number of Source Watersheds"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold")
  )
```

## Phosphorus Sensitivity Analysis

While phosphorus typically uses 100% efficiency, you can explore alternative scenarios:
```{r phosphorus_sensitivity, eval=FALSE}
# Standard phosphorus analysis
results_p_standard <- agri_classify_complete_custom(
  nugis_data,
  scale = "huc8",
  n_efficiency = 0.5,
  p_efficiency = 1.0
)

# Reduced phosphorus availability scenario
results_p_reduced <- agri_classify_complete_custom(
  nugis_data,
  scale = "huc8",
  n_efficiency = 0.5,
  p_efficiency = 0.8
)

# Compare phosphorus classifications
message("Standard (100% P efficiency):\n")
print(table(results_p_standard$P_class))

message("\nReduced Availability (80% P efficiency):\n")
print(table(results_p_reduced$P_class))
```

## Combined Sensitivity Analysis

Analyze both nutrients simultaneously with custom efficiencies:
```{r combined_sensitivity, eval=FALSE}
# Create multiple scenarios
scenarios <- list(
  baseline = list(n = 0.5, p = 1.0),
  improved_n = list(n = 0.7, p = 1.0),
  reduced_p = list(n = 0.5, p = 0.8),
  both_modified = list(n = 0.6, p = 0.9)
)

# Run all scenarios
results_list <- lapply(names(scenarios), function(scenario_name) {
  scenario <- scenarios[[scenario_name]]
  
  results <- agri_classify_complete_custom(
    nugis_data,
    scale = "huc8",
    n_efficiency = scenario$n,
    p_efficiency = scenario$p
  )
  
  # Return summary statistics
  data.frame(
    Scenario = scenario_name,
    N_Efficiency = scenario$n,
    P_Efficiency = scenario$p,
    N_Source = sum(results$N_class == "Source", na.rm = TRUE),
    P_Source = sum(results$P_class == "Source", na.rm = TRUE),
    N_Sink_Deficit = sum(results$N_class == "Sink_Deficit", na.rm = TRUE),
    P_Sink_Deficit = sum(results$P_class == "Sink_Deficit", na.rm = TRUE)
  )
})

# Combine results
sensitivity_summary <- do.call(rbind, results_list)
print(sensitivity_summary)
```

## Mapping Results Across Scenarios

Visualize spatial differences between efficiency scenarios:
```{r map_comparison, eval=FALSE}
# Join with spatial boundaries
spatial_standard <- boundaries %>%
  left_join(results_standard, by = c("huc8" = "ID"))

spatial_high <- boundaries %>%
  left_join(results_high, by = c("huc8" = "ID"))

# Create side-by-side maps
map_standard <- map_agricultural_classification(
  spatial_standard, 
  "nitrogen", 
  "N_class",
  "Standard Efficiency (50%)"
)

map_high <- map_agricultural_classification(
  spatial_high,
  "nitrogen",
  "N_class", 
  "Improved Efficiency (70%)"
)

# Display maps (requires cowplot or similar)
# cowplot::plot_grid(map_standard, map_high, ncol = 2)
```

## Interpreting Results

**Key Considerations:**

1. **Source Watersheds:** Higher efficiency factors typically increase the number of source watersheds as more nitrogen becomes "available" for redistribution

2. **Sink Classification:** Lower efficiencies may shift some "within watershed" areas to deficit status

3. **Management Implications:** Results inform decision-making about:
   - Required efficiency improvements to achieve redistribution goals
   - Sensitivity of regional planning to efficiency assumptions
   - Infrastructure needs under different management scenarios

4. **Uncertainty Quantification:** Range of results across scenarios indicates uncertainty due to efficiency assumptions

## Recommended Practice

For publication or policy analysis, consider:

1. **Report standard results** using established efficiency factors (0.5 for N, 1.0 for P)
2. **Conduct sensitivity analysis** with reasonable alternative values (e.g., ±20%)
3. **Document assumptions** clearly in methodology
4. **Discuss uncertainty** based on sensitivity results

## Example: Complete Sensitivity Workflow
```{r complete_workflow, eval=FALSE}
# Define efficiency range for nitrogen
n_efficiency_values <- seq(0.3, 0.7, by = 0.1)

# Run sensitivity analysis
sensitivity_results <- lapply(n_efficiency_values, function(eff) {
  results <- agri_classify_complete_custom(
    nugis_data,
    scale = "huc8",
    n_efficiency = eff,
    p_efficiency = 1.0
  )
  
  data.frame(
    N_Efficiency = eff,
    N_Source = sum(results$N_class == "Source", na.rm = TRUE),
    N_Sink_Deficit = sum(results$N_class == "Sink_Deficit", na.rm = TRUE),
    N_Within_Watershed = sum(results$N_class == "Within_Watershed", na.rm = TRUE)
  )
})

# Combine and visualize
sensitivity_df <- do.call(rbind, sensitivity_results)

# Create sensitivity plot
library(tidyr)
sensitivity_long <- sensitivity_df %>%
  pivot_longer(
    cols = c(N_Source, N_Sink_Deficit, N_Within_Watershed),
    names_to = "Classification",
    values_to = "Count"
  ) %>%
  mutate(
    Classification = gsub("N_", "", Classification),
    Classification = gsub("_", " ", Classification)
  )

ggplot(sensitivity_long, aes(x = N_Efficiency, y = Count, color = Classification)) +
  geom_line(size = 1.2) +
  geom_point(size = 3) +
  labs(
    title = "Nitrogen Classification Sensitivity to Efficiency Factor",
    x = "Nitrogen Efficiency Factor",
    y = "Number of Watersheds",
    color = "Classification"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    legend.position = "bottom"
  )
```

## Summary

Custom efficiency factors enable:
- **Sensitivity analysis** to quantify uncertainty
- **Scenario planning** for different management practices  
- **Regional adaptation** for varying agricultural systems
- **Policy evaluation** under different assumptions

For most applications, standard factors (N: 0.5, P: 1.0) remain appropriate, but sensitivity analysis provides valuable context for interpreting results and informing decision-making.

## References

- Spiegal et al. (2020). Manuresheds: Advancing nutrient recycling in US agriculture
- Flynn et al. (2023). Manureshed nutrient balances at watershed scales
- Akanbi et al. (2025). Towards circular nutrient economies: Integrated framework

For more information on standard manureshed methodologies, see `vignette("getting-started")`.