1 Introduction↩︎

Deep learning (DL) has become increasingly prevalent, and there is naturally a growing need for understanding DL predictions in many decision-making area, such as healthcare diagnosis and public policy-making. The high-stake nature of these areas means that these DL predictions are considered trustworthy only when they can be well explained. Meanwhile, time-series data has been frequently used in these areas [1]–[3], but it is always challenging to explain a time-series prediction due to the nature of temporal dependency and varying patterns over time. Moreover, time-series data often comes with confounding variables that affect both the input and output, making it even harder to explain predictions from DL models. On the other hand, many existing explanation methods are based on assigning importance scores for different parts of the input to explain model predictions  [4]–[9]. However, understanding the contribution of different input parts are usually not sufficiently informative for decision making: people often want to know what changes made to the input could have lead to a specific (desirable) prediction  [10]–[12]. We call such changed input that could have shifted the prediction to a specific target actionable counterfactual explanations. Below we provide an example in the context of time series.

Actionable counterfactual explanations help people understand how to achieve a counterfactual (target) output by modifying the current model input. However, such explanations may not be sufficiently informative in practice, especially under the causal effect of confounding variables which are often immutable. Specifically, some variables can hardly be changed once its value has been determined, and suggesting changing such variables are both meaningless and infeasible (e.g., a patient age and gender when modeling medical time series). This leads to a stronger requirement: a good explanation should make as few changes as possible on immutable variables; we call such explanations feasible counterfactual explanations. Below we provide an example in the context of time series.

In this paper, we propose a self-interpretable time series prediction model, dubbed Counterfactual Time Series (CounTS), which can both (1) perform time series predictions and (2) provide actionable and feasible counterfactual explanations for its predictions. Under common causal structure assumptions, our method is guaranteed to identify the causal effect between the input and output in the presence of exogenous (confounding) variables, thereby improving the generated counterfactual explanations’ feasibility. Our contribution is summarized as follows:

• We identify the actionability and feasibility requirements for generating counterfactual explanations for time series models and develop the first general self-interpretable method, dubbed CounTS, that satisfies such requirements.

• We provide theoretical guarantees that CounTS can identify the causal effect between the time series input and output in the presence of exogenous (confounding) variables, thereby improving feasibility in the generated explanations.

• Experiments on both synthetic and real-world datasets show that compared to state-of-the-art methods, CounTS significantly improves performance for generating counterfactual explanations while still maintaining comparable prediction accuracy.

2 Related Work↩︎

Interpretation Methods for Neural Networks. Various attribution-based interpretation methods have been proposed in recent years. Some methods focused on local interpretation [4]–[7], [13] while others are designed for global interpretation [14], [15]. The main idea is to assign attribution, or importance scores, to the input features in terms of their impact on the prediction (output). For example, such importance scores can be computed using gradients of the prediction with respect to the input [16]–[19]. Some interpretation methods are specialized for time series data; these include perturbation-based [20], rule-based [21], and attention-based methods [22], [23]. One typical method, Feature Importance in Time (FIT), evaluates the importance of the input data based on the temporal distribution shift and unexplained distribution shift [24]. However, these methods can only produce importance scores of the input features for the current prediction and therefore cannot generate counterfactual explanations (see Sec. 5 and Appendix 10 for empirical results).

Counterfactual Explanations for Time Series Models. There also works that generate counterfactual explanations for time series models. [25] proposed an association-rule algorithm to explain time series prediction by finding the frequent pairs of timestamps and generating counterfactual examples. [12] proposed a general explanation framework that generates counterfactual examples using residual generative adversarial networks (RGAN); it can be adapted for time series models. However, these works either fail to generate realistic counterfactual explanations (due to discretization error) or fail to generate feasible counterfactual explanations for time series models. In contrast, our CounTS as a principled variational causal method [26]–[28] can naturally generate realistic and feasible counterfactual explanations. Such advantages are empirically verified in Sec. 5.

Bayesian Deep Learning and Variational Autoencoders. Our work is also related to the broad categories of variational autoencoders (VAEs) [29] (which use inference networks to approximate posterior distributions) and Bayesian deep learning (BDL) [30]–[36] models (which use a deep component to process high-dimensional signals and a task-specific/graphical component to handle conditional/causal dependencies). [37] proposed the first VAE-based model for generating causal explanations for graph neural networks. [38], [39] proposed the first VAE-based models for performing causal inference and estimating treatment effect. However, none of them addressed the problem of counterfactual explanation, which involves solving an inverse problem to obtain the optimal counterfactual input. In contrast, our CounTS is the first VAE-based model to address this challenge, with theoretical guarantees and promising empirical results. From the perspective of BDL [31], [35], CounTS uses deep neural networks to process high-dimension signals (i.e., the deep component in [31]) and uses a Bayesian network to handle the conditional/causal dependencies among variables (i.e., the task-specific or graphical component in [31]). Therefore, CounTS is also the first BDL model for generating counterfactual explanations.

3 Preliminaries↩︎

Causal Model. Following the definition in [40], a causal model is described by a 3-tupple \(M = \left \langle U, V, F\right \rangle\). \(U\) is a set of exogenous variables \(\{u_1,\dots,u_m\}\) that is not determined by any other variables in this causal model. \(V\) is a set of endogenous variables \(\{v_1,\dots,v_n\}\) that are determined by variables in \(U \cup V\). We assume the causal model can be factorized according to a directed graph where each node represents one variable. \(F\) is a set of functions \(\{f_1,\dots,f_n\}\) describing the generative process of \(V\): \[v_i = f_i(pa_i, u_i), i=1,\dots,n,\] where \(pa_i\) denotes direct parent nodes of \(v_i\).

Counterfactual Inference. Counterfactual inference is interested in questions like “observing that \(X=x\) and \(Y=y\), what would be probability that \(Y=y^{cf}\) if the input \(X\) had been \(x^{cf}\)?". Formally, given a causal model \(\left \langle U, V, F\right \rangle\) where \(Y, X \in V\), counterfactual inference proceeds in three steps [40]:

Putting the three steps together, we have \[\begin{align} &P(Y_{X=x^{cf}}=y^{cf}|X=x,Y=y) \\ = &\sum\nolimits_u P(Y_{X=x^{cf}}(u)=y^{cf})P(u|x,y), \end{align}\] where \(x^{cf}\) can be an counterfactual explanation describing what would have shifted the outcome \(Y\) from \(y\) to \(y^{cf}\).

4 Method↩︎

In this section, we formalize the problem of counterfactual explanations for time series prediction and describe our proposed method for the problem.

Problem Setting. We focus on generating counterfactual explanations for predictions from time series models. We assume the model takes as input a multivariate time series \({\boldsymbol{x}}_i \in {\mathbb{R}}^{D \times T}\) and predicts the corresponding label \({\boldsymbol{y}}_i\), which can be a categorical label, a real value, or a time series \({\boldsymbol{y}}_i\in{\mathbb{R}}^{T}\). Given a specific input \({\boldsymbol{x}}_i\) and the model’s prediction \({\boldsymbol{y}}_i^{pred}\), our goal is to explain the model by finding a counterfactual time series \({\boldsymbol{x}}_i^{cf}\neq {\boldsymbol{x}}_i\) that could have lead the model to an alternative (counterfactual) prediction \({\boldsymbol{y}}_i^{cf}\).

4.1 Learning of CounTS↩︎

Causal Graph and Generative Model for CounTS. Our self-interpretable model CounTS is based on the causal graph in Fig. 1 (left), where \({\boldsymbol{x}}\in{\mathbb{R}}^{D\times T_{in}}\) is the input time series, \({\boldsymbol{y}}\) is the label, \({\boldsymbol{z}}\in{\mathbb{R}}^{H_z\times T_{mid}}\) is the representation of \({\boldsymbol{x}}\), \({\boldsymbol{u}}_l\in{\mathbb{R}}^{H_l\times T_{mid}}\) is the local exogenous variable, and \({\boldsymbol{u}}_g\in{\mathbb{R}}^{H_g}\) is the global exogenous variable. Both \({\boldsymbol{u}}_l\) and \({\boldsymbol{u}}_g\) are confounder variables; \({\boldsymbol{u}}_l\) can take different values in different time steps, while \({\boldsymbol{u}}_g\) is shared across all time steps.

This causal graph assumes the following factorization of the generative model \(p_\theta({\boldsymbol{y}},{\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})\): \[\begin{align} p_\theta({\boldsymbol{y}},{\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}}) = p_\theta({\boldsymbol{u}}_l,{\boldsymbol{u}}_g) p_\theta({\boldsymbol{z}}| {\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{x}})p_\theta({\boldsymbol{y}}| {\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{z}}),\label{eq:p95factor} \end{align}\tag{1}\] where \(\boldsymbol{\theta}\unboldmath\) denotes the collection of parameters for the generative model, and \(p_\theta({\boldsymbol{u}}_l,{\boldsymbol{u}}_g)=p_\theta({\boldsymbol{u}}_l)p_\theta({\boldsymbol{u}}_g)\). Here \(p_\theta({\boldsymbol{z}}| {\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{x}})\) and \(p_\theta({\boldsymbol{y}}| {\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{z}})\) are the encoder and the predictor, respectively. Specifically, we have \[\begin{align} p_\theta({\boldsymbol{u}}_l) &= {\mathcal{N}}({\boldsymbol{0}},{\boldsymbol{I}}),~~~~p_\theta({\boldsymbol{u}}_g) = {\mathcal{N}}({\boldsymbol{0}},{\boldsymbol{I}}), \\ p_\theta({\boldsymbol{z}}| {\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{x}}) &= {\mathcal{N}}(\mu_z({\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{x}}; \boldsymbol{\theta}\unboldmath), \sigma_z({\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{x}}; \boldsymbol{\theta}\unboldmath)), \\ p_\theta({\boldsymbol{y}}| {\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{z}}) &= {\mathcal{N}}(\mu_y({\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{z}}; \boldsymbol{\theta}\unboldmath), \sigma_y({\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}; \boldsymbol{\theta}\unboldmath)), \end{align}\] where \(\mu_z(\cdot; \cdot)\), \(\sigma_z(\cdot; \cdot)\), \(\mu_y(\cdot; \cdot)\), \(\sigma_y(\cdot; \cdot)\) are neural networks parameterized by \(\boldsymbol{\theta}\unboldmath\). Note that for classification models the predictor \(p_\theta({\boldsymbol{y}}| {\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{z}})\) becomes a categorical distribution \(Cat(f_y({\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{z}}; \boldsymbol{\theta}\unboldmath))\), where \(f_y(\cdot;\cdot)\) is a neural network.

Inference Model for CounTS. We use an inference model \(q_\phi({\boldsymbol{y}}, {\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})\) to approximate the posterior distribution of the latent variables, i.e., \(p_\theta({\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})\). As shown in Fig. 1 (right), we factorize \(q_\phi({\boldsymbol{y}}, {\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})\) as \[\begin{align} q_\phi({\boldsymbol{y}}, {\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}}) &= q_\phi({\boldsymbol{y}}|{\boldsymbol{x}})q_\phi({\boldsymbol{u}}_l,{\boldsymbol{u}}_g| {\boldsymbol{x}}, {\boldsymbol{y}})q_\phi({\boldsymbol{z}}| {\boldsymbol{x}}, {\boldsymbol{y}}),\label{eq:variational}\\ q_\phi({\boldsymbol{u}}_l,{\boldsymbol{u}}_g| {\boldsymbol{x}}, {\boldsymbol{y}})&=q_{\phi}({\boldsymbol{u}}_l|{\boldsymbol{x}}, {\boldsymbol{y}})q_{\phi}({\boldsymbol{u}}_g|{\boldsymbol{x}}, {\boldsymbol{y}}) \end{align}\tag{2}\] where \(\boldsymbol{\phi}\unboldmath\) is the collection of the inference model’s parameters, and . We parameterized each factor in Eqn. 2 as \[\begin{align} q_{\phi}({\boldsymbol{y}}|{\boldsymbol{x}}) &={\mathcal{N}}(\mu_y({\boldsymbol{x}}; \boldsymbol{\phi}\unboldmath), \sigma_y^2({\boldsymbol{x}}; \boldsymbol{\phi}\unboldmath)), \\ q_{\phi}({\boldsymbol{u}}_l|{\boldsymbol{x}}, {\boldsymbol{y}}) &= {\mathcal{N}}(\mu_{u_l}({\boldsymbol{x}}, {\boldsymbol{y}}; \boldsymbol{\phi}\unboldmath), \sigma_{u_l}^2({\boldsymbol{x}}, {\boldsymbol{y}}; \boldsymbol{\phi}\unboldmath)), \\ q_{\phi}({\boldsymbol{u}}_g|{\boldsymbol{x}}, {\boldsymbol{y}}) &= {\mathcal{N}}(\mu_{u_g}({\boldsymbol{x}}, {\boldsymbol{y}}; \boldsymbol{\phi}\unboldmath), \sigma_{u_g}^2({\boldsymbol{x}}, {\boldsymbol{y}}; \boldsymbol{\phi}\unboldmath)), \\ q_{\phi}({\boldsymbol{z}}|{\boldsymbol{x}}, {\boldsymbol{y}}) &= {\mathcal{N}}(\mu_z({\boldsymbol{x}}, {\boldsymbol{y}}; \boldsymbol{\phi}\unboldmath), \sigma_z^2({\boldsymbol{x}}, {\boldsymbol{y}}; \boldsymbol{\phi}\unboldmath)), \end{align}\] where \(\mu_{\cdot}(\cdot;\cdot)\) and \(\sigma_{\cdot}(\cdot;\cdot)\) denote neural networks with \(\boldsymbol{\phi}\unboldmath\) as their parameters.

Evidence Lower Bound. Our CounTS uses the evidence lower bound (ELBO) \({\mathcal{L}}_{ELBO}\) of the log likelihood \(\log p({\boldsymbol{y}}|{\boldsymbol{x}})\) as an objective to learn the generative and inference models. Maximizing the ELBO is equivalent to learning the optimal variational distribution \(q_\phi({\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})=\int q_\phi({\boldsymbol{y}}, {\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}}) d{\boldsymbol{y}}\) that best approximates the posterior distribution of the label and latent variables \(p_\theta({\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})\). Specifically, we have \[\begin{align} {\mathcal{L}}_{ELBO} = &~{\mathbb{E}}_{q_\phi({\boldsymbol{y}}, {\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})}[p_{\theta}({\boldsymbol{y}}, {\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})] \nonumber\\ &- {\mathbb{E}}_{q_\phi({\boldsymbol{y}}, {\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})}[q_\phi({\boldsymbol{y}}, {\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})].\label{eq:elbo95simple} \end{align}\tag{3}\] Note that different from typical ELBOs, we explicitly involves \({\boldsymbol{y}}\) and use \(q_\phi({\boldsymbol{y}}, {\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})\) rather than \(q_\phi({\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})\) (see Appendix 7 for details); this is to expose the factor \(q_\phi({\boldsymbol{y}}|{\boldsymbol{x}})\) in Eqn. 2 to allow for additional supervision on \({\boldsymbol{y}}\) (more details below). With the factorization in Eqn. 1 and Eqn. 2 , we can decompose Eqn. 3 as: \[\begin{align} {\mathcal{L}}&_{ELBO}=~{\mathbb{E}}_{q_\phi({\boldsymbol{y}}|{\boldsymbol{x}})}{\mathbb{E}}_{q_\phi({\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{z}}| {\boldsymbol{x}},{\boldsymbol{y}})}[p_\theta({\boldsymbol{y}}| {\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{z}})]\tag{4}\\ &-{\mathbb{E}}_{q_\phi({\boldsymbol{y}},{\boldsymbol{u}}_l,{\boldsymbol{u}}_g|{\boldsymbol{x}})}\big[KL[q_{\phi}({\boldsymbol{z}}|{\boldsymbol{x}},{\boldsymbol{y}})||p_\theta({\boldsymbol{z}}| {\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{x}})]\big]\tag{5}\\ &- {\mathbb{E}}_{q_\phi({\boldsymbol{y}}|{\boldsymbol{x}})}\big[KL[q_\phi({\boldsymbol{u}}_l,{\boldsymbol{u}}_g| {\boldsymbol{x}}, {\boldsymbol{y}})||p_\theta({\boldsymbol{u}}_l,{\boldsymbol{u}}_g)]\big]\tag{6}\\ &- {\mathbb{E}}_{q_{\phi}({\boldsymbol{y}}|{\boldsymbol{x}})}[q_{\phi}({\boldsymbol{y}}|{\boldsymbol{x}})],\tag{7} \end{align}\] where each term is computed with our neural network parameterization; Fig. 2 (left) shows the network structure. Below we briefly discuss the intuition of each term.

Final Objective Function. Inspired by [38], we include an additional term (\(N\) is the training set size) \[\begin{align} {\mathcal{L}}_y = \sum\nolimits_{i=1}^N\log q({\boldsymbol{y}}_i|{\boldsymbol{x}}_i) \end{align}\] to supervise \(q({\boldsymbol{y}}|{\boldsymbol{x}})\) using label \({\boldsymbol{y}}\) during training; this helps produce more accurate estimation for \({\boldsymbol{u}}_l\), \({\boldsymbol{u}}_g\), \({\boldsymbol{z}}\) using \(q_\phi({\boldsymbol{u}}_l,{\boldsymbol{u}}_g|{\boldsymbol{x}},{\boldsymbol{y}})\) and \(q_\phi({\boldsymbol{z}}|{\boldsymbol{x}},{\boldsymbol{y}})\). The final objective function then becomes \[\begin{align} {\mathcal{L}}_{CounTS} = {\mathcal{L}}_{ELBO} + \lambda {\mathcal{L}}_y,\label{eq:final} \end{align}\tag{8}\] where \(\lambda\) is a hyperparameter balancing both terms.

4.2 Inference Using CounTS↩︎

After learning both the generative model (Eqn. 1 ) and the inference model (Eqn. 2 ) by maximizing \({\mathcal{L}}_{CounTS}\) in Eqn. 8 , our self-interpretable CounTS can perform inference to (1) make predictions on \({\boldsymbol{y}}\) using the yellow branch in Fig. 2 (left), and (2) generate counterfactual explanation \({\boldsymbol{x}}^{cf}\) for any target label \({\boldsymbol{y}}^{cf}\) using the blue branch in Fig. 2 (right).

4.2.1 Prediction↩︎

We use the yellow branch in Fig. 2 (left) to predict \({\boldsymbol{y}}\). Specifically, given an input time series \({\boldsymbol{x}}\), the encoder will first infer the initial label \({\boldsymbol{y}}\) using \(q_\phi({\boldsymbol{y}}|{\boldsymbol{x}})\) and then infer \(({\boldsymbol{u}}_l,{\boldsymbol{u}}_g)\) and \({\boldsymbol{z}}\) from \({\boldsymbol{x}}\) and \({\boldsymbol{y}}\) using \(q_\phi({\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}},{\boldsymbol{y}})\). \(({\boldsymbol{u}}_l,{\boldsymbol{u}}_g)\) and \({\boldsymbol{z}}\) are then fed into the predictor \(p_\theta({\boldsymbol{y}}|{\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}})\) to infer the final label \({\boldsymbol{y}}\). Formally, CounTS predict the label as \[\begin{align} {\boldsymbol{y}}^{pred} = {\mathbb{E}}_{q_\phi({\boldsymbol{y}}|{\boldsymbol{x}})}{\mathbb{E}}_{q_\phi({\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{z}}| {\boldsymbol{x}},{\boldsymbol{y}})}[p_\theta({\boldsymbol{y}}| {\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{z}})]. \end{align}\] Empirically, we find that directly using the means of \(q_\phi({\boldsymbol{y}}|{\boldsymbol{x}})\) and \(q_\phi({\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{z}}| {\boldsymbol{x}},{\boldsymbol{y}})\) as input to \(p_\theta({\boldsymbol{y}}| {\boldsymbol{u}}_l,{\boldsymbol{u}}_g, {\boldsymbol{z}})\) already achieves satisfactory accuracy.

4.2.2 Generating Counterfactual Explanation↩︎

We use the blue branch in Fig. 2 (right) to generate counterfactual explanations via counterfactual inference. Our goal is to find the optimal counterfactual explanation \({\boldsymbol{x}}^{cf}\) defined below.

In words, we search for the optimal \({\boldsymbol{x}}^{cf}\) that would have shifted the model prediction from \({\boldsymbol{y}}^{pred}\) to \({\boldsymbol{y}}^{cf}\) while keeping \(({\boldsymbol{u}}_l,{\boldsymbol{u}}_g)\) unchanged. Since the definition of counterfactual explanations in Definition 1 involves causal inference with the intervention on \({\boldsymbol{x}}\), we need to first identify the causal probability \(p({\boldsymbol{y}}={\boldsymbol{y}}^{cf} | do({\boldsymbol{x}}={\boldsymbol{x}}^{\prime}), {\boldsymbol{u}})\) using observational probability, i.e., removing ‘do’ in the equation. The theorem below shows that this is achievable.

See Appendix 8.1 for the proof. With Theorem 1, we can rewrite Eqn. ?? as \[\begin{align} {\mathcal{L}}_{cf} = {\mathbb{E}}_{p({\boldsymbol{u}}| {\boldsymbol{x}}={\boldsymbol{x}}, {\boldsymbol{y}}={\boldsymbol{y}}^{pred})}{\mathbb{E}}_{p({\boldsymbol{z}}| {\boldsymbol{x}}^{\prime}, {\boldsymbol{u}})} [p({\boldsymbol{y}}^{cf} | {\boldsymbol{z}}, {\boldsymbol{u}}) ],\label{mqogbawd} \end{align}\tag{9}\] where \({\boldsymbol{u}}=({\boldsymbol{u}}_l,{\boldsymbol{u}}_g)\) and \(p({\boldsymbol{u}}| {\boldsymbol{x}}={\boldsymbol{x}}, {\boldsymbol{y}}={\boldsymbol{y}}^{pred})\) is approximated by \(q_\phi({\boldsymbol{u}}| {\boldsymbol{x}}={\boldsymbol{x}}, {\boldsymbol{y}}={\boldsymbol{y}}^{pred})\). We use Monte Carlo estimates to compute the expectation in Eqn. ?? and Eqn. 8 , iteratively compute the gradient \(\frac{\partial{\mathcal{L}}_{cf}}{\partial{\boldsymbol{x}}'}\) (via back-propagation) to search for the optimal \({\boldsymbol{x}}'\) in a way similar to [34], [41], and use it as \({\boldsymbol{x}}^{cf}\) (see the complete algorithm in Appendix 8.2).

5 Experiments↩︎

In this section, we evaluate our CounTS and existing methods on two synthetic and three real-world datasets. For each dataset, we evaluate different methods in terms of three metrics: (1) prediction accuracy, (2) counterfactual accuracy, and (3) counterfactual change ratio, with the last one as the most important metric. These metrics take different forms for different datasets (see details in Sec. 5.2-5.4).

5.1 Baselines and Implementations↩︎

We compare our CounTS with state-of-the-art methods for generating explanations for deep learning models, including Regularized Gradient Descent (RGD) [42], Gradient-weighted Class Activation Mapping (GradCAM) [16], Gradient SHapley Additive exPlanations (GradSHAP) [17], Local Interpretable Model-agnostic Explanations (LIME) [4], Feature Importance in Time (FIT) [24], Case-crossover APriori (CAP) [25], and Counterfactual Residual Generative Adversarial Network (CounteRGAN) [12] (see Appendix 9.2 for more details). Note that among these baselines, only RGD and CounteRGAN can generate actionable explanations. Other baselines, including FIT (which is designed for time series models), only provide importance scores as explanations; therefore some evaluation metrics may not be applicable for them (shown as ‘-’ in tables).

All methods above are implemented with PyTorch [43]. For fair comparison, the prediction model in all the baseline explanation methods has the same neural network architecture as the inference module in our CounTS. See the Appendix for more details on the architecture, training, and inference.

2pt

5.2 Toy Dataset↩︎

Dataset Description. We designed a toy dataset where the label is affected by only part of the input. A good counterfactual explanation should only modify this part of the input while keeping the other part unchanged. Following the causal graph in Fig. 1 (left), we have input \({\boldsymbol{x}}\in{\mathbb{R}}^{12}\) with each entry independently sampled from \(\mathcal{N} (\mu_x , \sigma^2_x)\) and an exogenous variable \(u \in {\mathbb{R}}\) (as the confounder) sampled from \(\mathcal{N} ( \mu_u , \sigma^2_u)\). To indicate which part of \({\boldsymbol{x}}\) affects the label, we introduce a mask vector \({\boldsymbol{m}}\in {\mathbb{R}}^{12}\) with first \(6\) entries \({\boldsymbol{m}}_{1:6}\) set to \(1\) and the last 6 entries \({\boldsymbol{m}}_{7:12}\) set to \(0\). We then generate \({\boldsymbol{z}}\in {\mathbb{R}}^{12}\) as \({\boldsymbol{z}}= {\boldsymbol{u}}\cdot({\boldsymbol{m}}\odot {\boldsymbol{x}})\) and the label \(y \sim Bern(\sigma({\boldsymbol{z}}^\top{\boldsymbol{1}} + u))\) where \(\sigma\) is the sigmoid function. Here \({\boldsymbol{x}}\)’s first \(6\) entries \({\boldsymbol{x}}_{1:6}\) is label-related and the last \(6\) entries \({\boldsymbol{x}}_{7:12}\) is label-agnostic.

Evaluation Metrics. We use three evaluation metrics:

Note that CCR is the most important metric among the three since our main focus is to generate actionable and feasible counterfactual explanations.

Quantitative Results. Table 1 compares our CounTS with the baselines in terms of the three metrics. CounTS outperforms all the baselines in terms of CCR with minimal prediction accuracy loss. This shows that our CounTS successfully identifies and fixes the exogenous variables \(u\) and \({\boldsymbol{m}}\) to generate actionable and feasible counterfactual explanations while still achieving prediction accuracy comparable to baselines. Note that actionable methods (i.e., CounTS, RGD, and CounteRGAN) outperforms importance score methods (i.e., FIT, LIME, and GradSHAP). This is expected because importance-score-based explanations can only explain the original prediction and therefore fail to infer what change on \({\boldsymbol{x}}\) will shift the prediction from \(y^{pred}\) to \(y^{cf}\). Our counterfactual accuracy is lower than CounteRGAN and RGD methods. This is reasonable since these baselines do not infer and fix the posterior distribution \(p({\boldsymbol{u}}|{\boldsymbol{x}},{\boldsymbol{y}})\) and are therefore more flexible for the generator (or back-propagation) to modify their input to push \(y^{pred}\) closer to the target \(y^{cf}\). However, such flexibility comes at the cost of low feasibility, reflected in their poor CCR performance.

2pt

5.3 Spike Dataset↩︎

Dataset Description. Inspired by [24], we construct the Spike synthetic dataset. In the dataset, each time series contains \(3\) channel, i.e., \({\boldsymbol{x}}\in{\mathbb{R}}^{3\times T}\). Each channel is a independent non–linear auto-regressive moving average (NARMA) sequence with randomly distributed spikes. The label sequence \({\boldsymbol{y}}\in{\mathbb{R}}^{T}\) starts at \(0\); it may switch to \(1\) as soon as there is a spike observed in any of the \(3\) channels and remain \(1\) until the last timestamp. A \(3\)-dimensional exogenous variable \({\boldsymbol{u}}\in\{0,1\}^{3}\) determines whether the spike in the channel can affect the final output; if so, we say the channel is active. Each of the three entries in \({\boldsymbol{u}}\), i.e., is sampled independently from three different Bernoulli distributions with parameters \(0.8\), \(0.4\), and \(0\), respectively (see more details on the dataset in Appendix 9.1).

Evaluation Metrics. We use three evaluation metrics:

Note that the scale of CCR will depend on the number of active channels, we therefore report our results for time series with \(1\) and \(2\) active channels. Higher CCR indicates better performance.

Quantitative Results. Table 2 shows the results for different methods in the Spike. Similar to the toy dataset, our CounTS outperforms all baselines in terms of CCR, and actionable methods (i.e., CounTS, RGD, and CounteRGAN) outperforms importance score methods (i.e., FIT, LIME, and GradSHAP) thanks to the former’s capability of modifying the input to shift the prediction towards the counterfactual target label.

Interestingly, besides promising performance in terms of CCR, our CounTS can also improve prediction performance, achieving lower prediction MSE. This is potentially due to CounTS’s ability to model the exogenous variable \({\boldsymbol{u}}\) that decides whether a spike in a specific channel affects the label \({\boldsymbol{y}}\). Similar to the toy datset, we notice both CounteRGAN and RGD methods have lower counterfactual MSE (both at \(0.074\)) than CounTS because they do not need to infer and fix the exogenous variable \({\boldsymbol{u}}\) (i.e., the mask) and therefore have more flexibility to modify the input \({\boldsymbol{x}}\). Since both CounteRGAN and RGD method are unaware of the mask, they suffer from lower CCR and produce more unwanted modification in inactive channels (more details below).

Qualitative Results. Fig. 3 shows an example time series as a case study. In this example, only channel \(1\) (blue) is active and thus the spike in channel \(1\) at timestamp \(t=16\) will flip the output \({\boldsymbol{y}}\) from \(0\) to \(1\). Both CounTS’s and RGD’s predictions \({\boldsymbol{y}}^{pred}\) (blue in Column 1) are very close to the ground truth.

Following the evaluation protocol above, we set our counterfactual (target) label \({\boldsymbol{y}}^{cf}\) by shifting \({\boldsymbol{y}}^{pred}\) by \(20\) timestamps to the right and padding \(0\)s on the left (yellow). We should expect an ideal counterfactual explanation to have no spike in \(t=[16, 36)\) and a new spike at \(t=36\) in channel 1. Since channel 2 and 3 are inactive, there should not be any changes. We can see all actionable methods (i.e., CounTS, CounteRGAN, and RGD) try to reduce the spikes in \(t=[16,36)\). However, CounteRGAN fails to remove the spike at \(t=29\); both RGD and CounteRGAN makes undesirable changes in inactive channels 2 and 3. Compared to baselines, our CounTS shows a significant advantage at inferring the active channel and suppressing the changes in inactive channels. This demonstrates that CounTS is more actionable and feasible than the existing time-series explanation methods. Note that FIT is not an actionable explanation method and therefore can only provide importance scores to explain the current prediction \({\boldsymbol{y}}^{pred}\) (see Appendix 10 Fig. 6 for details).

3.0pt

5.4 Real-World Datasets↩︎

Dataset Description. We evaluate our model on three real-world medical datasets: Sleep Heart Health Study (SHHS) [44], Multi-ethnic Study of Atherosclerosis (MESA) [45], and Study of Osteoporotic Fractures (SOF) [46], each containing subjects’ full-night breathing breathing signals. The breathing signals are divided into \(30\)-second segments; each segment has one of the sleep stages (‘Awake’, ‘Light Sleep’, ‘Deep Sleep’, ‘Rapid Eye Movement (REM)’) as its label. In total, there are \(2{,}651\), \(2{,}055\) and \(453\) patients in SHHS, MESA, and SOF, respectively, with an average of \(1{,}043\) breathing signal segments (approximately \(8.7\) hours).

Corresponding to our causal graph in Fig. 1 (left), \({\boldsymbol{x}}\) is the original breathing signal, \({\boldsymbol{z}}\) is signal patterns (representation) learned from the encoder, \({\boldsymbol{u}}\) can be ‘gender’, ‘age’, and even ‘dataset index’ (‘dataset index’ will be used during cross-dataset prediction; more details later), and \({\boldsymbol{y}}\) is the sleep stage label. Intuitively, ‘gender’, ‘age’, and ‘dataset index’ (due to different experiment instruments and signal measurement) can have causal effect on both the breathing signal patterns (subjects of different ages can have different breathing frequencies and magnitudes) and sleep stages (elderly subjects have less ‘Deep Sleep’ at night).

Evaluation Metrics. For brevity, we use abbreviations ‘A’, ‘R’, ‘L’, and ‘D’ to represent the four sleep stages ‘Awake’, ‘Rapid Eye Movement (REM)’, ‘Light Sleep’, ‘Deep Sleep’, respectively. We use three evaluation metrics:

In the experiments, we focus on cases where the target label \({\boldsymbol{y}}_{i,2}^{cf}\) is either ‘Awake’ or ‘Deep Sleep’ (two extremes of sleep stages), and evaluated \(6\) counterfactual action settings, namely A\(\rightarrow\)D, L\(\rightarrow\)D, R\(\rightarrow\)D, R\(\rightarrow\)A, D\(\rightarrow\)A, and L\(\rightarrow\)A.

Cross-Dataset and Intra-Dataset Settings. Since the three datasets are collected with different devices and procedures, we consider the dataset index as an additional exogenous variable in \({\boldsymbol{u}}\), which has causal effect on both the learned representation \({\boldsymbol{z}}\) and the sleep stage \({\boldsymbol{y}}\). This leads to two different settings: the intra-dataset setting (without the dataset index confounder) and the cross-dataset setting (with the dataset index confounder).

Quantitative Results. Table 3 shows the average CCR of different methods for both intra-dataset and cross-dataset settings. Results show that our CounTS outperforms the baselines in all the dataset settings in terms of the average CCR. This demonstrates CounTS’s capability of generating feasible and actionable explanations in complex real-world datasets. Table 4 shows the detailed CCR for each counterfactual action setting in SHHS, where CounTS is leading in almost all counterfactual action settings. Detailed results for other datasets can be found in Table 7\(\sim\)11 of Appendix 10.

2pt

1pt

Table 5 shows the prediction accuracy for both intra-dataset and cross-dataset settings. Similar to Table 1 and Table 2, all baselines (e.g., RGD, CounteRGAN, and FIT) explain the same prediction model and therefore share the same prediction accuracy. Results show that CounTS achieves prediction accuracy comparable to the original prediction model.  Table 6 shows the average counterfactual accuracy for different dataset settings over \(6\) counterfactual action settings. As in the toy and Spike datasets, RGD achieves higher counterfactual accuracy because it does not need to infer and fix the exogenous variable \({\boldsymbol{u}}\) and hence enjoys more flexibility when modifying the input \({\boldsymbol{x}}\); this often leads to undesirable changes (e.g., breathing frequency and amplitude which are related to the subject’s age), less feasible counterfactual explanations, and therefore worse CCR performance. In contrast, breathing frequency and amplitude are potentially captured by \({\boldsymbol{u}}_l\) and \({\boldsymbol{u}}_g\) in CounTS and kept unchanged (see qualitative results below). Detailed counterfactual accuracy for every counterfactual action setting and every dataset setting can be found in Appendix 10.

Qualitative Analysis: Problem Setting. As some background on sleep staging, note that a patient’s breathing signal will be much more periodic in ‘Deep Sleep’ than in ‘Awake’. Fig. 4 shows an example time series from the MESA dataset as a case study. It contains two consecutive segments of breathing signals (\(60\) seconds in total). The first segment (\(0\sim 30\) seconds) is regular and periodic, and therefore both CounTS and baselines correctly predict its label as ‘Deep Sleep’ (‘D’); the second segment (\(30\sim 60\) seconds) is irregular and therefore both CounTS and baselines correctly predict its label as ‘Awake’ (‘A’). The goal is to generate a counterfactual explanation such that the model predicts the second segment as ‘Deep Sleep’.

Qualitative Analysis: Ideal Counterfactual Explanation. Since a ‘Deep Sleep’ segment should be more periodic, an ideal counterfactual explanation should keep the first segment (\(0\sim 30\) seconds) unchanged and maintain its periodicity, make the second segment (\(30\sim 60\) seconds) more periodic (the pattern for ‘Deep Sleep’), and keep the breathing frequency (number of breathing cycles) unchanged throughout both segments (this is related to "feasibility" since a patient usually has the same breathing frequency during ‘Deep Sleep’).

Qualitative Analysis: Detailed Results. Fig. 4 compares the generated counterfactual explanations from CounTS, RGD, and CounteRGAN, and we can see that our CounTS’s explanation is closer to the ideal case. Specifically:

6 Conclusion↩︎

In this paper, we identified the actionability and feasibility requirements for time series models counterfactual explanations and proposed the first self-interpretable time series prediction model, CounTS, that satisfies both requirements. Our theoretical analysis shows that CounTS is guaranteed to identify the causal effect between time series input and output in the presence of confounding variables, thereby generating counterfactual explanations in the causal sense. Empirical results showed that our method achieves competitive prediction accuracy on time series data and is capable of generating more actionable and feasible counterfactual explanations. Interesting future work includes further reducing the computation complexity during counterfactual inference, handling uncertainty in the explanation [47], and extending our method to multimodal data.

Acknowledgement↩︎

We would like to thank the reviewers/AC for the constructive comments to improve the paper. JY and HW are partially supported by NSF Grant IIS-2127918. We also thank National Sleep Research Resource (NSRR) data repository and the researchers for sharing SHHS [48], [49], MESA [48], [50] and SOF [48], [51].

langley00?

7 Evidence Lower Bound↩︎

We can derive the evidence lower bound (ELBO) in Eqn. 3 as follows: \[\begin{align} &\log p({\boldsymbol{y}}|{\boldsymbol{x}}) \\ &=\log \int_{{\boldsymbol{z}}}\int_{{\boldsymbol{u}}} p({\boldsymbol{y}}, {\boldsymbol{z}},{\boldsymbol{u}}| {\boldsymbol{x}}) d {\boldsymbol{u}}d {\boldsymbol{z}}\\ &\geq \mathbb{E}_{q_\phi ({\boldsymbol{u}}, {\boldsymbol{z}}, {\boldsymbol{y}}| {\boldsymbol{x}})} \left [ \log\frac{p_\theta({\boldsymbol{y}}, {\boldsymbol{u}}, {\boldsymbol{z}}| {\boldsymbol{x}})}{q_\phi({\boldsymbol{u}}, {\boldsymbol{z}}, {\boldsymbol{y}}| {\boldsymbol{x}})}\right ] \\ &=\mathbb{E}_{q_\phi({\boldsymbol{u}}, {\boldsymbol{z}}, {\boldsymbol{y}}| {\boldsymbol{x}})} \left [\log\frac{p({\boldsymbol{u}}) p_\theta({\boldsymbol{z}}| {\boldsymbol{u}}, {\boldsymbol{x}})p_\theta({\boldsymbol{y}}| {\boldsymbol{u}}, {\boldsymbol{z}})}{q_\phi({\boldsymbol{u}}, {\boldsymbol{z}}, {\boldsymbol{y}}| {\boldsymbol{x}})} \right ] \\ &={\mathbb{E}}_{q_\phi({\boldsymbol{y}}, {\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})}[p_{\theta}({\boldsymbol{y}}, {\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})] - {\mathbb{E}}_{q_\phi({\boldsymbol{y}}, {\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})}[q_\phi({\boldsymbol{y}}, {\boldsymbol{u}}_l,{\boldsymbol{u}}_g,{\boldsymbol{z}}|{\boldsymbol{x}})]\\ &=\mathbb{E}_{q_\phi({\boldsymbol{u}}, {\boldsymbol{z}}, {\boldsymbol{y}}| {\boldsymbol{x}})}[\log p({\boldsymbol{u}})] +\mathbb{E}_{q_\phi({\boldsymbol{u}}, {\boldsymbol{z}}, {\boldsymbol{y}}| {\boldsymbol{x}})}[\log p_\theta({\boldsymbol{z}}| {\boldsymbol{u}}, {\boldsymbol{x}})] +\mathbb{E}_{q_\phi({\boldsymbol{u}}, {\boldsymbol{z}}, {\boldsymbol{y}}| {\boldsymbol{x}})}[\log p_\theta({\boldsymbol{y}}| {\boldsymbol{u}}, {\boldsymbol{z}})] -\mathbb{E}_{q_\phi({\boldsymbol{u}}, {\boldsymbol{z}}, {\boldsymbol{y}}| {\boldsymbol{x}})}[\log q_\phi({\boldsymbol{u}}, {\boldsymbol{z}}, {\boldsymbol{y}}| {\boldsymbol{x}})] \label{eq:marginal} \end{align}\tag{10}\]

8 More Details on Counterfactual Inference↩︎

8.1 Proof on Identifiability↩︎

8.2 Counterfactual Explanation Algorithm↩︎

The pseudo-code for counterfactual explanation is shown in Algorithm 5.

9 More Details on Experiments↩︎

9.1 Details on Datasets↩︎

Spike Dataset. The generation process of the Spike dataset is summarized below:

We generate \(D=3\) independent channels of non–linear auto-regressive moving average (NARMA) time series data using the following formula: \[x_{d,t+1}=0.5x_{d,t}+0.5x_{d,t}\sum_{i=0}^{l-1}x(t-l)+1.5u(t-(l-1))u(t)+0.5+\alpha_{d}t \label{eq:spike95narma}\tag{11}\] for \(t=[1,\dots,80]\), order \(l=2\), \(u \sim {\mathcal{N}}(0, 0.03)\), and \(\alpha_d\) is set differently for each channel (\(\alpha_1 = 0.1\), \(\alpha_2 = 0.065\) and \(\alpha_3 = 0.003\)). We use \(d\) to index the \(3\) channels.

\[\begin{align} &\mathbf{\boldsymbol{\theta}\unboldmath} = [0.8, 0.4, 0] ; \\ &n_{d} \sim \operatorname{Bernoulli}(\mathbf{\theta}_{d}) ; \\ &m_{d} \sim \operatorname{Bernoulli}(\mathbf{\theta}_{d}) ; \\ &\eta_{d}= \begin{cases}\operatorname{Poisson}(\lambda=2) & \mathbf{1} \text{ if } \left(n_{d}==1\right) \\ 0 & \text{ otherwise }\end{cases} \\ &\mathbf{g}_{d} = \operatorname{Sample}\left([T], \eta_{d}\right)\\ &x_{d, t}=x_{d, t}+(\kappa_{d,t}+ \theta_{d}) \quad \text{where } \kappa_{d,t}\sim {\mathcal{N}}(1, 0.3) \quad \forall t \in \mathbf{g}_{d} \label{synthetic:2} \\ & y_{t}=\left\{\begin{array}{ll} 0 & t \leqslant \min \left(\mathbf{g}_d\right) \text{where } m_d=1 \\ 1 & \text{otherwise} \end{array}\right. \end{align}\tag{12}\]

Real-World Datasets. The real-world datasets also include additional patient information such as age, gender and race. The age range is \([44, 90]\) for SHHS, \([54, 95]\) for MESA and \([71, 90]\) for SOF. In total, there are \(2{,}651\), \(2{,}055\) and \(453\) patients in SHHS, MESA, and SOF, respectively, with an average of \(1{,}043\) breathing signal segments (approximately \(8.7\) hours).

9.2 Details on Baselines↩︎

We use the following five types of state-of-the-art baselines:

• Gradient-based methods. Regularized Gradient Descent (RGD) [42] directly models \(p(y|x)\) and provide the explanation by modifying input with gradients along with L1 regularization; it is therefore it is an actionable explanation method. Gradient-weighted Class Activation Mapping (GradCAM) [16] is originally designed to models with convolutional layers. We added convolutional layers to our model to adapt GradCAM for our time series data. Gradient SHapley Additive exPlanations (GradSHAP) [17] is a game theoretic approach that uses the expectation of gradients to approximate the SHAP values.

• Perturbation-based method. Local Interpretable Model-agnostic Explanations (LIME) [4] is a local interpretation approach based on the local linearity assumption and provides explanations by fitting the output of the model given locally perturbed input.

• Distribution shift method. Feature Importance in Time (FIT) [24] is a time-series back-box model explanation method that evaluates the importance of the input data based on the temporal distribution shift and unexplained distribution shift.

• Association rule method. Case-crossover APriori (CAP) [25] applies association rule mining algorithm, Apriori, to explore the causal relationship in time-series data.

• Generative method. Counterfactual Residual Generative Adversarial Network (CounteRGAN) [12] combines a Residual GAN (RGAN) and a classifier to generate counterfactual output. The generated residual output is considered as the result of the \(do\) operation. CounteRGAN is an actionable explanation method.

10 Additional Results↩︎

Results on the Spike Dataset. Note that methods such as FIT are not an actionable explanation method. It cannot provide counterfactual explanation \({\boldsymbol{x}}^{cf}\) that could shift the model prediction from \({\boldsymbol{y}}^{pred}\) to \({\boldsymbol{y}}^{cf}\); it can only provide importance scores to explain the current prediction \({\boldsymbol{y}}^{pred}\). Fig. 6 shows the importance scores produced by FIT to explain the same model in Fig. 3 in the Spike dataset.

Results on Real-World Datasets. As FIT cannot provide actionable explanations, the importance scores produced by FIT to explain the same model in Fig. 4 in real-world dataset are separately shown in Fig. 7.

Table 7\(\sim\)11 show more detailed CCR results with different counterfactual action settings and dataset settings. Table ¿tbl:tab:realcfaccu? shows detailed counterfactual accuracy for each dataset settings and each counterfactual action settings.

References↩︎