.blue-dark[.mediumlarge[Characterizing climate pathways using echo state networks and feature importance]]

.title[
# .blue-dark[.mediumlarge[Characterizing climate pathways using echo state networks and feature importance]]
]
.author[
### Katherine Goode
]
.author[
### April 3, 2024
]
.author[
### Montana State University Statistics Department Seminar
]
.author[
### .dark-blue[.smaller[SAND2024-04036PE]]
]

---

.pull-left-small{
  width: 35%;
  float: left;
}
.pull-right-large{
  width: 60%;
  float: right;
}
.pull-right-small{
  width: 35%;
  float: right;
}
.pull-left-large{
  width: 65%;
  float: left;
}
.pull-left-smallish{
  width: 46%;
  float: left;
}
.pull-right-largeish{
  width: 53%;
  float: right;
}
.pull-left-largeish{
  width: 53%;
  float: left;
}
.pull-right-smallish{
  width: 46%;
  float: right;
}
.pull-right-smaller{
  width: 20%;
  float: right;
}
.pull-left-larger{
  width: 76%;
  float: left;
}
.pull-left-smaller{
  width: 22%;
  float: left;
}
.pull-right-larger{
  width: 76%;
  float: right;
}
.pull-left-v3 {
  width: 60%;
  height: 92%;
  float: left;
}
.pull-right-v3 {
  width: 35%;
  float: right;
  padding-left: 1%;
}

</style>

## But first... Working as a statistician at a national lab

### .blue[Sandia National Laboratories]

Federally Funded Research and Development Center (FFRDC)

From the [Sandia website](https://www.sandia.gov/):

> "We strive to become the laboratory that the U.S. turns to first for technology solutions to the most challenging problems that threaten peace and freedom for our nation and the globe."

]

Recent labs news:

<img src="figs/lab-news.png" width="100%" style="display: block; margin: auto;" />
]

---

## But first... Working as a statistician at a national lab

### .blue[Sandia National Laboratories]

Federally Funded Research and Development Center (FFRDC)

From the [Sandia website](https://www.sandia.gov/):

> "We strive to become the laboratory that the U.S. turns to first for technology solutions to the most challenging problems that threaten peace and freedom for our nation and the globe."

]

Examples of programs at Sandia:

<img src="figs/labs-programs.png" width="90%" style="display: block; margin: auto;" />
]

---

## Previous Projects

**Predicting Explosive Device Characteristics**

]

]

---

## Previous Projects

**Determining Processing Conditions of Nuclear Particulates**

---

## Current Projects

**Different Perspectives on Trust/Explainability/Interpretability with Machine Learning**

.pull-left-large[
- Random forest explainability
  - Identify tree topology patterns in random forest
  
- Private yet explainable
  - Develop differentially private + explainable clustering algorithm

- Maturity levels for SciML
  - Develop framework for assessing maturity of machine learning models used to complement computational simulation models
  
- Trust in AI
  - Human study to quantify trust in an algorithm when transparency and interactivity of algorithm are adjusted
]

]

---

# Back to Climate Project

---

## Remainder of Talk...

- .blue[Overview]: Project Motivation and Goals

- .blue[Approach]: Echo State Networks and Feature Importance

- .blue[Climate Application]: Mount Pinatubo

- .blue[Ending]: Conclusions and Future/Current Work

Acknowledgements

- **Sub-thrust team**: Daniel Ries, Kellie McClernon

- **Thrust lead**: Lyndsay Shand

- **General advise**: Gabriel Huerta, J. Derek Tucker

---

# Overview

### .bright-teal[Project Motivation and Goals]

---

## Climate Interventions

.center[.smallmedium[Image source: [https://eos.org/science-updates/improving-models-for-solar-climate-intervention-research](https://eos.org/science-updates/improving-models-for-solar-climate-intervention-research)]]

]

- Stratospheric aerosol injections

- Marine cloud brightening

- Cirrus cloud thinning

- etc.
]

---

## CLDERA Grand Challenge

---

## Climate Event Exemplar

- Mount Pinatubo eruption in 1991

- Released 18-19 Tg of sulfur dioxide

- Proxy for anthropogenic stratospheric aerosol injection

.pull-left[
<img src="figs/pinatubo-ash-cloud.jpg" width="90%" style="display: block; margin: auto;" />
.center[.smallmedium[Image source: [https://pubs.usgs.gov/fs/1997/fs113-97/](https://pubs.usgs.gov/fs/1997/fs113-97/)]]
]

.pull-right[
 
<img src="figs/pinatubo.jpg" width="90%" style="display: block; margin: auto;" />
 
.center[.smallmedium[Image source: [https://volcano.si.edu/volcano.cfm?vn=273083](https://volcano.si.edu/volcano.cfm?vn=273083)]]
]

---

## Observational Thrust

**Objective**: Develop algorithms to .blue[characterize (i.e., quantify) relationships between climate variables] related to a climate event using observational data

.pull-left-smaller[
<img src="figs/obs-data-collection.png" width="60%" style="display: block; margin: auto;" />
]

.pull-right-larger[
<img src="figs/obs-pathway.png" width="100%" style="display: block; margin: auto;" />
 
.center[.smallmedium[Image credit: CLDERA leadership]]
]

---

## Mount Pintabuo Pathway

AOD increased as a result of injection of sulfur dioxide [1; 2]]

.medium[Temperatures at pressure levels of 30-50 mb rose 2.5-3.5 degrees centigrade compared to 20-year mean [3]] 
]

---

## Mount Pintabuo Pathway

.center[.smaller[Figure generated using Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA- 2) data [4].]]

---

## Our Approach

Use machine learning...

**Step 1: Model climate event variables with echo state network**

Allow complex machine learning model to capture complex variable relationships

<img src="figs/esn-step1.png" width="95%" style="display: block; margin: auto;" />
]

**Step 2: Quantify relationships via explainability**

Apply feature importance to understand relationships captured by model

<img src="figs/esn-step2.png" width="100%" style="display: block; margin: auto;" />
]

---

# Approach

### .bright-teal[Step 1: Echo State Networks]

---

## Echo State Networks

**Overview**

- Nonlinear machine learning model for temporal data
- Sibling to recurrent neural network (RNN)
- Computationally efficient due to many parameters being randomly sampled instead of estimated

.blue[Output stage]: ridge regression
`$$\textbf{y}_{t} = \mathbf{V} \mathbf{h}_t + \boldsymbol{\epsilon}_{t} \ \ \ \ \ \ {\bf \epsilon_t } \sim N(\textbf{0}, \sigma^2_\epsilon \textbf{I})$$`
.blue[Hidden stage]: nonlinear stochastic transformation
`$$\mathbf{h}_t = g_h \left(\frac{\nu}{|\lambda_w|} \mathbf{W} \mathbf{h}_{t-1} + \mathbf{U} \mathbf{\tilde{x}}_{t-\tau}\right)$$`
`$$\tilde{\mathbf{x}}_{t-\tau}=\left[\textbf{x}'_{t-\tau},\textbf{x}'_{t-\tau-\tau^*},...,\mathbf{x}'_{t-\tau-m\tau^*}\right]'$$`
]

.pull-right[
 
 
 
Elements of `$\textbf{W}$` and `$\textbf{U}$` randomly sampled as:
`\begin{align}
 \textbf{W}[h,c_w] &=\gamma_{h,c_w}^w\mbox{Unif}(-a_w,a_w)+(1-\gamma_{h,c_w}^w)\delta_0\\
 \textbf{U}[h,c_u] &=\gamma_{h,c_u}^u\mbox{Unif}(-a_u,a_u)+(1-\gamma_{h,c_u}^u)\delta_0
\end{align}`
where `$\gamma_{h,c_w}^w \sim Bern(\pi_w)$`, `$\gamma_{h,c_u}^u \sim Bern(\pi_u)$`, `$\delta_0$` is a Dirac function, and values of `$a_w$`, `$a_u$`, `$\pi_w$`, and `$\pi_u$` are pre-specified and set to small values.
]

---

## Echo State Networks

---

## Echo State Networks: Spatio-Temporal Context

Recall that we are working with spatio-temporal data...

---

## Echo State Networks: Spatio-Temporal Context

Spatio-temporal processes at spatial locations `$\{\textbf{s}_i\in\mathcal{D}\subset\mathbb{R}^2;i=1,...,N\}$` over times `$t=1,...,T$`...

.pull-left[
.blue[Output variable] (e.g., stratospheric temperature): 
  
`$${\bf Z}_{Y,t} = \left(Z_{Y,t}({\bf s}_1),Z_{Y,t}({\bf s}_2),...,Z_{Y,t}({\bf s}_N)\right)'$$`

]

.blue[Input variables] (e.g., lagged aerosol optical depth and stratospheric temperature): For `$k=1,...,K$`
  
`$${\bf Z}_{k,t} = \left(Z_{k,t}({\bf s}_1),Z_{k,t}({\bf s}_2),...,Z_{k,t}({\bf s}_N)\right)'$$`
]

| Stage | Formula | Description |
| ----- | ------- | ----------- |
| Output data stage | `${\bf Z}_{Y,t}\approx\boldsymbol{\Phi}_Y\textbf{y}_{t}$` | Basis function decomposition (e.g., PCA) |
| Output stage | `$\textbf{y}_{t} = \mathbf{V} \mathbf{h}_t + \boldsymbol{\epsilon}_{t}$` | Ridge regression |
| Hidden stage | `$\mathbf{h}_t = g_h \left(\frac{\nu}{\lvert\lambda_w\rvert} \mathbf{W} \mathbf{h}_{t-1} + \mathbf{U} \mathbf{\tilde{x}}_{t-\tau}\right)$` `$\tilde{\mathbf{x}}_{t-\tau}=\left[\textbf{x}'_{t-\tau},\textbf{x}'_{t-\tau-\tau^*},...,\mathbf{x}'_{t-\tau-m\tau^*}\right]'$` | Nonlinear stochastic transformation |
| Input data stage | `${\bf Z}_{k,t}\approx\boldsymbol{\Phi}_k\textbf{x}_{k,t}$` `$\ \ \ \ \ \ \textbf{x}_t=[\textbf{x}'_{1,t},...,\textbf{x}'_{K,t}]'$` | Basis function decomposition (e.g., PCA) |

---

## Echo State Networks: Spatio-Temporal Context

---

# Approach

### .bright-teal[Step 2: Feature Importance]

---

## Feature Importance for ESNs

**Goal:** Quantify the effect of .blue[input variable] `$k$` over .blue[block of times] `$(t-b,...,t-1)$` on .blue[forecasts] at time `$t$`

---

## Feature Importance for ESNs

.pull-right-small[
**Concept**: "Adjust" inputs at times(s) of interest and quantify effect on model performance

- .blue[Permute values]: spatio-temporal permutation feature importance (stPFI)

- .blue[Set values to zero]: spatio-temporal zeroed feature importance (stZFI)
]

**Feature Importance**: Difference in RMSEs from "adjusted" and observed spatial predictions:

`$$RMSE_{adj,t}-RMSE_{obs,t}$$`

**Interpretation**: Large feature importance indicates "adjusted" inputs lead to a decrease in model performance indicating the model finds those inputs important for prediction

---

## Feature Importance for ESNs

Computing feature importance of variable 1 at time `$t=5$` using a block size of 3:

`$$RMSE_{adj,5}-RMSE_{obs,5}$$`

---

# Climate Application

### .bright-teal[Mount Pinatubo]

---

## Mount Pinatubo Example: Data

**Data** Monthly values from 1980 to 1995 and -86 to 86 degrees latitude

**Source** Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA- 2)

`$$\frac{Z_{k,mth,yr}({\bf s}_i)-\bar{Z}_{k,mth,\cdot}({\bf s}_i)}{sd(Z_{k,mth,\cdot}({\bf s}_i))}$$`

where `$\bar{Z}_{k,mth,\cdot}({\bf s}_i)$` and `$sd(Z_{k,mth,\cdot}({\bf s}_i))$` are average and standard deviation at location `${\bf s}_i$` during month `$mth$` for variable `$k$`
]

.pull-right-larger[
<img src="figs/merra2_weighted_global_averages.png" width="90%" style="display: block; margin: auto;" />
]

---

## Mount Pinatubo Example: Model

.pull-left-small[
**Goal** Forecast stratospheric temperature (50mb) one month ahead given lagged values of stratospheric temperature (50mb) and AOD

**Details**

- All variables: First 5 principal components computed on normalized anomalies

- Hyperparameters selected via grid search when years of 1994 and 1995 were held out

- Fit 25 ESNs to account for random variability
]

<img src="figs/merra2_rmse.png" width="97%" style="display: block; margin: auto;" />
]

---

## Mount Pinatubo Example: Feature Importance

Block size of 6 months (including variability across 25 ESNs)

---

## Mount Pinatubo Example: Feature Importance

Observed and adjusted RMSEs with block size of 6 months (including variability across 25 ESNs)

---

## Mount Pinatubo Example: Feature Importance

Block sizes of 1-6 months (averaged over 25 ESNs)

---

# Conclusions and Future Work

---

## Summary and Conclusions

- Interested in quantifying relationships between climate variables associated with pathway of climate event

- Motivated by increasing possibility of climate interventions

- Our machine learning approach:

- Use ESN to model variable relationships

- Understand variable relationships using proposed spatio-temporal feature importance

- Approach provided evidence of AOD being a variable in Mount Pinatubo climate pathway affecting stratospheric temperature

---

## Future (Current) Work

**Methodology**

- Ensembles of ESNs

- Implement proposed retraining technique [5] to lessen detection of spurious relationships due to correlation

- Feature importance computed on spatial scale

- Comparison to other newly proposed explainability techniques for ESNs (layer-wise relevance propagation)  [6]

**Mt. Pinatubo Application**

- Inclusion of additional Mt. Pinatubo pathway variables (e.g., SO2, radiative flux, surface temperature)

- Importance of grouped variables

- Application to climate simulation ensembles

- Use of climate simulation ensembles for method assessment

---

## References

.smallmedium[
[1] S. Guo, G. J. Bluth, W. I. Rose, et al. "Re-evaluation of SO$_2$
release of the 15 June 1991 Pinatubo eruption using ultraviolet and
infrared satellite sensors". In: _Geochemistry, Geophysics, Geosystems_
5 (4 2004), pp. 1-31. DOI:
[10.1029/2003GC000654](https://doi.org/10.1029%2F2003GC000654).

[2] M. Sato, J. E. Hansen, M. P. McCormick, et al. "Stratospheric
aerosol optical depths, 1850-1990". In: _Journal of Geophysical
Research: Atmospheres_ 98.D12 (1993), pp. 22987-22994. DOI:
[https://doi.org/10.1029/93JD02553](https://doi.org/https%3A%2F%2Fdoi.org%2F10.1029%2F93JD02553).
eprint:
https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/93JD02553. URL:
[https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/93JD02553](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/93JD02553).

[3] K. Labitzke and M. McCormick. "Stratospheric temperature increases
due to Pinatubo aerosols". In: _Geophysical Research Letters_ 19 (2
1992), pp. 207-210. DOI:
[10.1029/91GL02940](https://doi.org/10.1029%2F91GL02940).

[4] R. Gelaro, W. McCarty, M. J. Suarez, et al. "The ModernEra
Retrospective Analysis for Research and Applications, Version 2
(MERRA-2)". In: _Journal of Climate_ 30 (14 2017), pp. 5419-5454. DOI:
[10.1175/JCLI-D-16-0758.1](https://doi.org/10.1175%2FJCLI-D-16-0758.1).

[5] G. Hooker, L. Mentch, and S. Zhou. "Unrestricted permutation forces
extrapolation: variable importance requires at least one more model, or
there is no free variable importance". In: _Statistics and Computing_
31 (2021), pp. 1-16.

[6] M. Landt-Hayen, P. Kröger, M. Claus, et al. "Layer-Wise Relevance
Propagation for Echo State Networks Applied to Earth System
Variability". In: _Signal, Image Processing and Embedded Systems
Trends_. Ed. by D. C. Wyld. Computer Science & Information Technology
(CS & IT): Conference Proceedings 20. ARRAY(0x55588c8d8680), 2022, pp.
115-130. ISBN: 978-1-925953-80-0. DOI:
[doi:10.5121/csit.2022.122008](https://doi.org/doi%3A10.5121%2Fcsit.2022.122008).
URL:
[https://doi.org/10.5121/csit.2022.122008](https://doi.org/10.5121/csit.2022.122008).
]

---

# Thank you

.white[Slides: [goodekat.github.io/posts/2024-04-03.html](https://goodekat.github.io/posts/2024-04-03.html)]

---

# Back-Up Slides

---

## Weighted RMSE

Weighted RMSE to compute feature importance (weighted by cos latitude):

`$$\mbox{weighted RMSE}_t =\sqrt{\frac{\sum_{loc}w_{loc}(y_{t,loc}-\hat{y}_{t,loc,adj})^2}{\sum_{loc}w_{loc}}} - \sqrt{\frac{\sum_{loc}w_{loc}(y_{t,loc}-\hat{y}_{t,loc})^2}{\sum_{loc}w_{loc}}}$$`

---

## Feature Importance: Spatio-Temporal Context

Contribution of spatial locations to ZFI: `$\sqrt{\frac{w_{loc}(y_{t,loc}-\hat{y}_{t,loc,zeroed})^2}{\sum_{loc}w_{loc}}} - \sqrt{\frac{w_{loc}(y_{t,loc}-\hat{y}_{t,loc})^2}{\sum_{loc}w_{loc}}}$`

---

## Feature Importance: Spatio-Temporal Context

Contribution of spatial locations to ZFI: `$\sqrt{\frac{w_{loc}(y_{t,loc}-\hat{y}_{t,loc,zeroed})^2}{\sum_{loc}w_{loc}}} - \sqrt{\frac{w_{loc}(y_{t,loc}-\hat{y}_{t,loc})^2}{\sum_{loc}w_{loc}}}$`

---

## Simulated Data Demonstration

**Simulated response**

`$$Z_{Y,t}({\bf s}_i)=Z_{2,t}({\bf s}_i) \beta + \delta_t({\bf s}_i) + \epsilon_t({\bf s}_i)$$`

where

- `$Z_{2,t}$` spatio-temporal covariate
- `$\delta_t({\bf s}_i)$` spatio-temporal random effect
- `$\epsilon_t({\bf s}_i) \overset{iid}{\sim}  N(0,\sigma_{\epsilon}^2)$`

]

<img src="figs/syn_data_heat.png" width="90%" style="display: block; margin: auto;" />
]

---

## Simulated Data: Effect of Variability on FI

---

## Effect of Correlation on FI