---
title: "xaiHydro: Explainable AI Tools for Hydro-Climate Modelling"
subtitle: "A Unified R Workflow — Book Chapter Companion Vignette"
author:
  - name: "Sadikul Islam"
    affiliation: "Department of Hydrology, [Your Institution], India"
    email: "sadikul.islam@institution.ac.in"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 4
    number_sections: true
    fig_width: 8
    fig_height: 5
    df_print: kable
vignette: >
  %\VignetteIndexEntry{xaiHydro: Explainable AI for Hydro-Climate Models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: references.bib
link-citations: true
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.align = "center",
  warning   = FALSE,
  message   = FALSE,
  fig.width  = 8,
  fig.height = 5,
  out.width  = "100%"
)
```

---

> **Citation**: Islam, S. (2026). Explainable AI for Hydro-Climate Models:
> A Unified R Workflow Using SHAP, LIME and Partial Dependence. In:
> *Hydro-Climate Analytics: Remote Sensing, AI and Geospatial Modelling*
> (Springer). DOI: 10.XXXX/XXXXXXX

---

# Introduction

## Background

Machine learning models are increasingly used in hydro-climate science
to predict streamflow, estimate drought indices, map flood risk, and
project evapotranspiration under changing climate scenarios. These
models — random forests, gradient boosting machines, deep neural networks
— often outperform traditional physically-based models in predictive
accuracy [@breiman2001random]. However, they operate as **black boxes**:
their internal decision logic is opaque to practitioners, regulators, and
policy makers who must trust model outputs for operational water
management and climate adaptation decisions.

Explainable Artificial Intelligence (XAI) addresses this gap by providing
mathematically rigorous tools to attribute model predictions to individual
input features. XAI transforms a black-box model into an interpretable
system without sacrificing predictive performance.

## The xaiHydro package

`xaiHydro` provides a **unified, hydrology-aware R workflow** for applying
three complementary XAI families to any trained hydro-climate model:

| Method | Scope | Key advantage |
|--------|-------|---------------|
| **SHAP** [@lundberg2017unified] | Global + local | Theoretically grounded; satisfies efficiency, symmetry, and dummy axioms |
| **LIME** [@ribeiro2016should] | Local | Interpretable surrogate; works with any model type |
| **PDP / ALE** [@friedman2001greedy; @apley2020visualizing] | Global marginal | Visualises feature response curves across the predictor range |

All three methods require only a **trained model** and its **training
data** — no additional field data collection, no satellite downloads, no
pre-processing pipelines. This makes `xaiHydro` uniquely accessible for
researchers who already have a fitted model and want to explain it.

## Scope of this vignette

This vignette serves as both:

1. **A package tutorial** — covering every function in `xaiHydro` with
   reproducible code.
2. **The book chapter case study** — the outputs produced here are the
   figures and analyses referenced in the accompanying Springer chapter.

---

# Installation

```r
# CRAN release (recommended)
install.packages("xaiHydro")

# Development version from GitHub
# install.packages("remotes")
remotes::install_github("sadikul-islam/xaiHydro")
```

Load the package:

```{r load-pkg}
library(xaiHydro)
```

---

# Simulated datasets

`xaiHydro` ships with two built-in data simulators so every example in
this vignette is fully reproducible without external data.

## Daily streamflow dataset

```{r sim-streamflow}
df_flow <- sim_streamflow_data(n = 365, seed = 2026)
head(df_flow, 6)
```

```{r sim-str}
str(df_flow)
```

The function generates a synthetic daily hydro-meteorological dataset
mimicking realistic seasonal dynamics:

| Column | Description | Unit |
|--------|-------------|------|
| `precipitation` | Daily rainfall | mm |
| `temperature` | Mean daily air temperature | °C |
| `evapotranspiration` | Reference ET (Penman–Monteith style) | mm/day |
| `soil_moisture` | Volumetric soil moisture | m³/m³ |
| `snow_water_equiv` | Snow water equivalent | mm |
| `ndvi` | NDVI vegetation greenness index | dimensionless |
| `lag1_precip` | Previous-day precipitation | mm |
| `lag3_precip` | 3-day moving mean precipitation | mm |
| `streamflow` | **Response variable** — simulated discharge | m³/s |

The response is generated as a non-linear function of predictors with
added noise, including a precipitation × NDVI interaction term to
capture vegetation-modulated runoff.

## Monthly drought index dataset

```{r sim-drought-show}
df_drought <- sim_drought_data(n = 240, seed = 2026)
head(df_drought, 6)
```

This dataset simulates 20 years of monthly hydro-climatic anomalies with
SPEI (Standardised Precipitation-Evapotranspiration Index) as the
response, driven by precipitation anomalies, temperature anomalies, SPI
at 3- and 12-month scales, PDSI, NDVI anomaly, and ET deficit.

---

# Model training

## Random Forest — streamflow prediction

We train a Random Forest model on the simulated streamflow dataset.
The same workflow applies to any regression model supported by `predict()`.

```{r train-rf, eval=FALSE}
library(randomForest)

# Separate predictors and response
y_flow <- df_flow$streamflow
X_flow <- df_flow[, setdiff(names(df_flow), "streamflow")]

set.seed(42)
rf_model <- randomForest(
  x         = X_flow,
  y         = y_flow,
  ntree     = 300,
  mtry      = 3,
  importance = TRUE
)

print(rf_model)
# Mean of squared residuals: ~X.XX
# % Var explained:            ~XX.X
```

## XGBoost — drought index prediction

```{r train-xgb, eval=FALSE}
library(xgboost)

y_spei   <- df_drought$spei
X_drought <- df_drought[, setdiff(names(df_drought), "spei")]

xgb_mat   <- xgboost::xgb.DMatrix(as.matrix(X_drought), label = y_spei)
xgb_model <- xgboost::xgboost(
  data      = xgb_mat,
  nrounds   = 150,
  max_depth = 4,
  eta       = 0.1,
  objective = "reg:squarederror",
  verbose   = 0
)
```

---

# The hydro explainer

## Creating an explainer

`hydro_explainer()` is the single entry point for all XAI analyses.
It wraps [DALEX::explain()] [@biecek2021explanatory] and attaches
hydrology-specific metadata used for axis labels, titles, and
interpretation messages throughout the package.

```{r explainer, eval=FALSE}
exp_flow <- hydro_explainer(
  model    = rf_model,
  data     = X_flow,
  y        = y_flow,
  variable = "streamflow",
  units    = "m3/s",
  label    = "RandomForest"
)

print(exp_flow)
# == xaiHydro Explainer ==
#   Target variable : streamflow [m3/s]
#   Model label     : RandomForest
#   Observations    : 365
#   Predictors      : 8
#   Model class     : randomForest
```

The explainer object is compatible with all `xaiHydro` functions and
can also be passed directly to any `DALEX` function for additional
analyses.

## Explainer for XGBoost (custom predict wrapper)

For models that require a non-standard prediction call (XGBoost, Keras,
PyTorch via `torch`), supply a `predict_function`:

```{r explainer-xgb, eval=FALSE}
pfun_xgb <- function(model, newdata) {
  predict(model, xgboost::xgb.DMatrix(as.matrix(newdata)))
}

exp_drought <- hydro_explainer(
  model            = xgb_model,
  data             = X_drought,
  y                = y_spei,
  variable         = "drought_index",
  units            = "SPEI",
  label            = "XGBoost",
  predict_function = pfun_xgb
)
```

---

# SHAP Analysis

SHAP (SHapley Additive exPlanations) values decompose each model
prediction into additive contributions from each feature, satisfying
three desirable axioms: efficiency (contributions sum to the prediction
gap from baseline), symmetry (equal features receive equal credit), and
dummy (non-contributing features get zero) [@lundberg2017unified].

## Computing SHAP values

```{r shap-compute, eval=FALSE}
shap_vals <- hydro_shap(
  explainer = exp_flow,
  nsim      = 50,     # Monte Carlo replicates — increase for stability
  seed      = 42
)

# Returns a tidy data frame:
# id | feature | shap_value | feature_value | mean_abs_shap
head(shap_vals, 8)
```

**Tip:** `nsim = 50` is adequate for publication figures. Use
`nsim = 100–200` for final manuscript outputs.

## Global SHAP summary plot

```{r shap-summary, eval=FALSE}
plot_shap_summary(shap_vals, top_n = 8)
```

The beeswarm summary plot shows:

- **Y-axis** — features ranked by mean |SHAP| (most important at top).
- **X-axis** — SHAP value (positive = pushes prediction up; negative =
  pulls it down).
- **Colour** — raw feature value (red = high, blue = low).

Features where red points cluster on the right indicate that high
values of that predictor increase predicted streamflow — physically
consistent with rainfall–runoff theory.

## Single-observation waterfall

```{r shap-waterfall, eval=FALSE}
# Explain a high-flow day (e.g., observation 180 after a storm event)
plot_shap_waterfall(shap_vals, obs_id = 180, top_n = 7)
```

The waterfall plot shows exactly how prediction for day 180 was built
up from the global baseline (mean predicted streamflow) through additive
feature contributions. Red bars push the prediction above baseline;
blue bars pull it below.

---

# LIME Explanation

LIME (Local Interpretable Model-agnostic Explanations) fits a locally
weighted linear surrogate model in the neighbourhood of a specific
observation [@ribeiro2016should]. It provides a simple, interpretable
approximation of the black-box model's local behaviour.

```{r lime, eval=FALSE}
lime_result <- hydro_lime(
  explainer    = exp_flow,
  new_obs      = X_flow[180, , drop = FALSE],
  n_features   = 6,      # top-6 features in the local model
  kernel_width = 0.75
)

# Local R² indicates surrogate model quality (closer to 1 = better)
plot_lime_hydro(lime_result)
```

The local R² is printed in the subtitle. Values above 0.7 indicate the
surrogate captures the model's local structure well.

---

# Partial Dependence Profiles

Partial Dependence Plots (PDP) show the **marginal effect** of a single
predictor on the response, averaging over all other predictors
[@friedman2001greedy]. Accumulated Local Effects (ALE) are preferred
when predictors are correlated, as they avoid averaging over unrealistic
covariate combinations [@apley2020visualizing].

## PDP — marginal effects

```{r pdp, eval=FALSE}
pdp_result <- hydro_pdp(
  explainer   = exp_flow,
  variable    = c("precipitation", "soil_moisture",
                   "lag1_precip",  "temperature"),
  type        = "partial",
  grid_points = 60
)

plot_pdp_hydro(pdp_result, ncol = 2)
```

## ALE — preferred for correlated predictors

```{r ale, eval=FALSE}
ale_result <- hydro_pdp(
  explainer = exp_flow,
  variable  = c("precipitation", "soil_moisture"),
  type      = "accumulated"
)

plot_pdp_hydro(ale_result, ncol = 2)
```

In hydrology, precipitation and soil moisture are typically positively
correlated (wet soils follow wet periods). ALE correctly handles this
by computing effects within observed data slices rather than averaging
over impossible combinations.

---

# Variable Importance

Permutation importance measures how much model error increases when
each feature's values are randomly shuffled, breaking the feature–
response relationship [@breiman2001random]. Features causing large
error increases are the most important.

```{r importance, eval=FALSE}
imp <- hydro_importance(
  explainer     = exp_flow,
  loss_function = "rmse",
  B             = 20       # repetitions — reduces variance of estimate
)

plot_importance_hydro(imp, top_n = 8)
```

---

# Prediction Breakdown

The breakdown profile decomposes a single prediction into ordered
additive contributions from each feature [@biecek2021explanatory].
Unlike SHAP (which averages over all orderings), the sequential Break
Down algorithm shows contributions in a specific order; the SHAP-type
breakdown averages over `B` random orderings for a more stable result.

```{r breakdown, eval=FALSE}
bd <- hydro_breakdown(
  explainer = exp_flow,
  new_obs   = X_flow[180, , drop = FALSE],
  type      = "shap",   # or "break_down"
  B         = 25
)

plot_breakdown_hydro(bd, max_features = 8)
```

---

# Residual Diagnostics

Before trusting XAI explanations, model residuals should be examined.
Systematic patterns in residuals (heteroscedasticity, temporal
autocorrelation) can bias SHAP and PDP estimates.

```{r residuals, eval=FALSE}
hydro_residuals(exp_flow)
```

The four-panel figure checks:

1. **Residual histogram** — departure from normality (skewness, heavy
   tails).
2. **Residuals vs predicted** — heteroscedasticity; systematic bias at
   high flows.
3. **Observed vs predicted** — overall fit; 1:1 line is perfect
   prediction.
4. **Residuals over index** — temporal autocorrelation in daily data
   (important for hydrological models with lag effects).

---

# Full automated XAI report

`hydro_xai_report()` runs all key analyses in one call and assembles
them into a publication-ready multi-panel figure. This is the figure
used in the book chapter case study.

```{r report, eval=FALSE}
hydro_xai_report(
  explainer = exp_flow,
  nsim      = 50,
  top_n     = 6,
  new_obs   = X_flow[180, , drop = FALSE],
  seed      = 42,
  save_path = "xaiHydro_streamflow_report.png",
  width     = 16,
  height    = 11,
  dpi       = 300
)
```

---

# Complete drought index case study

This section demonstrates the full `xaiHydro` workflow on an XGBoost
drought model, reproducing the second case study in the book chapter.

```{r drought-full, eval=FALSE}
library(xgboost)

# Data
df_d  <- sim_drought_data(n = 240, seed = 2026)
y_d   <- df_d$spei
X_d   <- df_d[, setdiff(names(df_d), "spei")]

# Model
xgb_m <- xgboost::xgboost(
  data      = xgboost::xgb.DMatrix(as.matrix(X_d), label = y_d),
  nrounds   = 150, max_depth = 4, eta = 0.1,
  objective = "reg:squarederror", verbose = 0
)

# Explainer
pfun  <- function(m, nd) predict(m, xgboost::xgb.DMatrix(as.matrix(nd)))
exp_d <- hydro_explainer(xgb_m, X_d, y_d,
                          variable         = "drought_index",
                          units            = "SPEI",
                          label            = "XGBoost",
                          predict_function = pfun)

# XAI suite
shap_d <- hydro_shap(exp_d, nsim = 50)
imp_d  <- hydro_importance(exp_d, B = 20)
pdp_d  <- hydro_pdp(exp_d, variable = c("spi_3", "spi_12", "et_deficit"))

plot_shap_summary(shap_d, top_n = 6,
                  title = "SHAP Summary — SPEI drought model (XGBoost)")
plot_importance_hydro(imp_d)
plot_pdp_hydro(pdp_d, ncol = 3)
hydro_residuals(exp_d)
```

---

# Interpreting XAI results in a hydro-climate context

## Physical interpretation of SHAP

When `precipitation` and `lag1_precip` appear in the top-ranked SHAP
features with positive effects (red points on the right), this confirms
the model has correctly learnt that higher antecedent and current
rainfall drives higher streamflow — consistent with rainfall–runoff
theory. Conversely, high `evapotranspiration` having negative SHAP
(blue points on left) indicates that high ET demand reduces net water
available for runoff.

## ALE vs PDP selection guide

Use **ALE** when:
- Predictors are correlated (precipitation–soil moisture, temperature–ET).
- The training data does not uniformly cover the predictor space.
- You suspect extrapolation artefacts in PDP.

Use **PDP** when:
- Predictors are approximately independent.
- You want a simpler figure for a non-technical audience.

## LIME local R² threshold

An R² < 0.5 for the LIME surrogate means the local linear approximation
is unreliable for that observation. Increase `kernel_width` (to use a
broader neighbourhood) or reduce `n_features` (for a sparser model).

---

# Session information

```{r session}
sessionInfo()
```

---

# References