Inside the LSTM: An XAI Field Guide to Weather Prediction

LSTMs are still the go-to architecture for a lot of time series work, but they're annoying to trust. You get a number out the other end and no real sense of why the model landed there. This tutorial walks through training an LSTM on daily temperature data, then pulling it apart with three explainability methods: permutation importance, SHAP, and Integrated Gradients.

Who this is for: people who already know some Keras and want to add interpretability to a forecasting model, not a from-scratch intro to neural nets.

LSTMs want a 3D tensor — (samples, timesteps, features)

— so before anything else we need to turn a flat column of temperatures into overlapping 7-day windows, each one paired with the value on day 8.

import numpy as np
import pandas as pd
from sklearn.preprocessing import MinMaxScaler

df = pd.read_csv("weather_data.csv")
data = df['Temperature'].values.reshape(-1, 1)

scaler = MinMaxScaler(feature_range=(0, 1))
scaled_data = scaler.fit_transform(data)

X, y = [], []
for i in range(7, len(scaled_data)):
    X.append(scaled_data[i-7:i])
    y.append(scaled_data[i])

X, y = np.array(X), np.array(y)
print(f"Input shape: {X.shape}") # Output: (Samples, 7, 1)

Scaling matters more than it sounds like it should — LSTMs trained on unscaled temperature values are prone to exploding gradients, and training just falls apart. The windowing step is really the whole trick here: every prediction only ever sees the past seven days, nothing more.

Two stacked LSTM layers, dropout after each one, early stopping so we don't have to babysit the epoch count.

from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Dropout, Input
from tensorflow.keras.callbacks import EarlyStopping

model = Sequential([
    Input(shape=(7, 1)),
    LSTM(units=100, activation='relu', return_sequences=True),
    Dropout(0.2),
    LSTM(units=100, activation='relu'),
    Dropout(0.2),
    Dense(units=1)
])

model.compile(optimizer='adam', loss='mse')

early_stopping = EarlyStopping(
    monitor='val_loss', 
    patience=10, 
    restore_best_weights=True
)

history = model.fit(
    X, y, 
    epochs=50, 
    batch_size=32, 
    validation_split=0.2, 
    callbacks=[early_stopping]
)

One detail worth flagging: return_sequences=True

on the first LSTM layer is required, not optional, if you're stacking a second one — it passes the full sequence of hidden states forward instead of just the final one. patience=10

gives the model some room to wander before we give up on it, and restore_best_weights=True

means we walk away with the checkpoint from its best epoch, not whatever it happened to end on.

Once training is done, the obvious question is: does the model lean on yesterday's temperature, or does it weight the whole week roughly evenly? Three different methods, three different angles on that question.

The idea is blunt but effective — shuffle one lag day across the whole dataset, see how much worse the predictions get, repeat for each of the 7 days.

from sklearn.metrics import mean_squared_error

base_preds = model.predict(X, verbose=0)
base_error = mean_squared_error(y, base_preds)

feature_importance = []
for i in range(7):
    X_permuted = X.copy()
    np.random.shuffle(X_permuted[:, i, 0]) # Shuffle a specific lag day

    permuted_preds = model.predict(X_permuted, verbose=0)
    permuted_error = mean_squared_error(y, permuted_preds)

    feature_importance.append(permuted_error - base_error)

Whichever day causes the biggest jump in error when scrambled is the one the model depends on most. In practice, for a weather series, that's usually the most recent day or two — but it's worth checking rather than assuming.

Permutation importance tells you what matters on average across the dataset. SHAP tells you what mattered for one specific prediction, which is a lot more useful when you're trying to explain a single forecast to someone. KernelExplainer

works reasonably well for LSTMs, though be warned it's slow — it's model-agnostic and treats the network as a black box, so it has to make a lot of prediction calls.

import shap

background_data = X[:100].reshape(-1, 7)
instance_to_explain = X[0:1].reshape(-1, 7)

def model_predict_wrapper(input_2d):
    return model.predict(input_2d.reshape(-1, 7, 1), verbose=0)

explainer = shap.KernelExplainer(model_predict_wrapper, background_data)
shap_values = explainer.shap_values(instance_to_explain)

shap.initjs()
shap.force_plot(
    explainer.expected_value[0], 
    shap_values[0].flatten(), 
    feature_names=[f'Day {i+1}' for i in range(7)],
    features=instance_to_explain[0,:]
)

The force plot shows, day by day, whether each lag pushed that particular forecast above or below the baseline average. It's the kind of chart you can actually put in front of a non-technical stakeholder.

This one's gradient-based rather than perturbation-based, and it plays nicely with TensorFlow since it just needs access to the gradient tape. You interpolate a path from a baseline (usually all zeros) to the real input, and integrate the gradients along that path.

import tensorflow as tf

def integrated_gradients(model, input_tensor, baseline=None, steps=50):
    if baseline is None:
        baseline = tf.zeros_like(input_tensor)

    alphas = tf.linspace(0.0, 1.0, steps + 1)
    interpolated_inputs = [baseline + alpha * (input_tensor - baseline) for alpha in alphas]
    interpolated_inputs = tf.concat(interpolated_inputs, axis=0)

    @tf.function
    def call_model(inputs):
        return model(inputs, training=False)

    with tf.GradientTape() as tape:
        tape.watch(interpolated_inputs)
        predictions = call_model(interpolated_inputs)

    grads = tape.gradient(predictions, interpolated_inputs)

    avg_grads = (grads[:-1] + grads[1:]) / 2.0
    return tf.reduce_mean(avg_grads, axis=0) * (input_tensor - baseline)

sample = tf.convert_to_tensor(X[0:1], dtype=tf.float32)
ig_attributions = integrated_gradients(model, sample)

The steps

parameter is a tradeoff — more steps gives a smoother, more accurate integral approximation, at the cost of more forward passes. 50 is a reasonable default; you rarely need to go much higher for a 7-day window.

None of these three methods is strictly "correct" — they answer slightly different questions. Permutation importance gives you the global picture, SHAP gives you a per-prediction breakdown that's easy to explain to someone else, and Integrated Gradients gives you a gradient-based view that's cheap to compute since it doesn't need to touch the model as a black box. Running all three and comparing where they agree (and where they don't) tends to be more informative than picking just one.