# Self-Interpretable Time Series Prediction with Counterfactual Explanations

Jingquan Yan  Hao Wang

###### Abstract

Interpretable time series prediction is crucial for safety-critical areas such
as healthcare and autonomous driving. Most existing methods focus on
interpreting predictions by assigning important scores to segments of time
series. In this paper, we take a different and more challenging route and aim
at developing a self-interpretable model, dubbed Counterfactual Time Series
(CounTS), which generates counterfactual and actionable explanations for time
series predictions. Specifically, we formalize the problem of time series
counterfactual explanations, establish associated evaluation protocols, and
propose a variational Bayesian deep learning model equipped with
counterfactual inference capability of time series abduction, action, and
prediction. Compared with state-of-the-art baselines, our self-interpretable
model can generate better counterfactual explanations while maintaining
comparable prediction accuracy. Code will be available at
https://github.com/Wang-ML-Lab/self-interpretable-time-series.

Machine Learning, ICML

##  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 (Zhao et al., 2021; Jin et al., 2022; Yang et
al., 2022), 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 (Ribeiro et al., 2016; Lundberg & Lee, 2017b; Chen et al., 2018;
Wang et al., 2019b; Weinberger et al., 2020; Plumb et al., 2020). 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
(Wachter et al., 2017b; Goyal et al., 2019; Nemirovsky et al., 2022). 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.

######  Example 1 (Actionable Counterfactual Explanation).

Suppose there is a model that takes as input a time series of breathing signal
${\bf x}\in{\mathbb{R}}^{T}$ from a subject of age $u=60$ to predict the
corresponding sleep stage as
$y^{pred}=\mbox{`Awake'}\in\\{\mbox{`Awake'},\mbox{`Light Sleep'},\mbox{`Deep
Sleep'}\\}$. Typical methods assign importance scores to each entry of ${\bf
x}$ to explain the prediction. However, they do not provide _actionable_
counterfactual explanations on how to modify ${\bf x}$ to ${\bf x}^{cf}$ such
that the prediction can change to $y^{cf}=\mbox{`Deep Sleep'}$. An ideal
method with such capability could provide more information on why the model
make specific predictions.

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.

######  Example 2 (Feasible Counterfactual Explanation).

In Example 1, age $u$ is a confounder that affects both ${\bf x}$ and $y$
since elderly people (i.e., larger $u$) are more likely to have irregular
breathing ${\bf x}$ and more ‘Awake’ time (i.e., $y=\mbox{`Awake'}$) at night.
To generate a counterfactual explanation to change $y^{pred}$ to ‘Deep Sleep’,
typical methods tend to suggest decreasing the age $u$ from $60$ to $50$,
which is _infeasible_ (since age cannot be changed in practice). An ideal
method would first infer the age $u$ and search for a _feasible_
counterfactual explanation ${\bf x}^{cf}$ that could change $y^{pred}$ to
‘Deep Sleep’ while keeping $u$ unchanged.

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 (Ribeiro et al., 2016; Lundberg & Lee, 2017b;
Plumb et al., 2018; Chen et al., 2018; Wang et al., 2019b) while others are
designed for global interpretation (Ghorbani et al., 2019; Natesan Ramamurthy
et al., 2020). 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 (Selvaraju et al., 2017; Lundberg & Lee,
2017a; Shrikumar et al., 2017; Sundararajan et al., 2017). Some interpretation
methods are specialized for time series data; these include perturbation-based
(Pan et al., 2021), rule-based (Rajapaksha & Bergmeir, 2022), and attention-
based methods (Heo et al., 2018; Lim et al., 2021). 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
(Tonekaboni et al., 2020). 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 D for empirical
results).

Counterfactual Explanations for Time Series Models. There also works that
generate counterfactual explanations for time series models. (Dhaou et al.,
2021) proposed an association-rule algorithm to explain time series prediction
by finding the frequent pairs of timestamps and generating counterfactual
examples. (Nemirovsky et al., 2022) 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 (Wang et al., 2020; Mao et al., 2021b; Gupta et al., 2021) 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) (Kingma & Welling,
2013) (which use inference networks to approximate posterior distributions)
and Bayesian deep learning (BDL) (Wang et al., 2015; Wang & Yeung, 2016; Wang,
2017; Huang et al., 2019; Wang et al., 2019a; Wang & Yeung, 2020; Ding et al.,
2022) models (which use a deep component to process high-dimensional signals
and a task-specific/graphical component to handle conditional/causal
dependencies). (Lin et al., 2022) proposed the first VAE-based model for
generating causal explanations for graph neural networks. (Louizos et al.,
2017; Pawlowski et al., 2020) 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 (Wang & Yeung, 2016, 2020), CounTS uses deep neural
networks to process high-dimension signals (i.e., the deep component in (Wang
& Yeung, 2016)) and uses a Bayesian network to handle the conditional/causal
dependencies among variables (i.e., the task-specific or graphical component
in (Wang & Yeung, 2016)). Therefore, CounTS is also the first BDL model for
generating counterfactual explanations.

##  3 Preliminaries

Causal Model. Following the definition in (Pearl, 2009), 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 (Pearl, 2009):

  1. 1.

Abduction. Calculate the posterior distribution of $u$ given the observation
$X=x$ and $Y=y$, i.e., $P(u|X=x,Y=y)$.

  2. 2.

Action. Perform causal intervention on the variable $X$, i.e., $do(X=x^{cf})$.

  3. 3.

Prediction. Calculate the counterfactual probability
$P(Y_{X=x^{cf}}(u)=y^{cf})$ with respect to the posterior distribution of
$P(u|X=x,Y=y)$.

Putting the three steps together, we have

$\displaystyle P(Y_{X=x^{cf}}=y^{cf}|X=x,Y=y)$ $\displaystyle=$
$\displaystyle\sum\nolimits_{u}P(Y_{X=x^{cf}}(u)=y^{cf})P(u|x,y),$

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 ${\bf x}_{i}\in{\mathbb{R}}^{D\times T}$ and predicts
the corresponding label ${\bf y}_{i}$, which can be a categorical label, a
real value, or a time series ${\bf y}_{i}\in{\mathbb{R}}^{T}$. Given a
specific input ${\bf x}_{i}$ and the model’s prediction ${\bf y}_{i}^{pred}$,
our goal is to explain the model by finding a counterfactual time series ${\bf
x}_{i}^{cf}\neq{\bf x}_{i}$ that could have lead the model to an alternative
(counterfactual) prediction ${\bf y}_{i}^{cf}$.

Figure 1: Left: The causal graph and generative model for CounTS. Right:
Inference model for CounTS.

###  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 ${\bf
x}\in{\mathbb{R}}^{D\times T_{in}}$ is the input time series, ${\bf y}$ is the
label, ${\bf z}\in{\mathbb{R}}^{H_{z}\times T_{mid}}$ is the representation of
${\bf x}$, ${\bf u}_{l}\in{\mathbb{R}}^{H_{l}\times T_{mid}}$ is the local
exogenous variable, and ${\bf u}_{g}\in{\mathbb{R}}^{H_{g}}$ is the global
exogenous variable. Both ${\bf u}_{l}$ and ${\bf u}_{g}$ are confounder
variables; ${\bf u}_{l}$ can take different values in different time steps,
while ${\bf u}_{g}$ is shared across all time steps.

This causal graph assumes the following factorization of the generative model
$p_{\theta}({\bf y},{\bf u}_{l},{\bf u}_{g},{\bf z}|{\bf x})$:

$\displaystyle p_{\theta}({\bf y},{\bf u}_{l},{\bf u}_{g},{\bf z}|{\bf
x})=p_{\theta}({\bf u}_{l},{\bf u}_{g})p_{\theta}({\bf z}|{\bf u}_{l},{\bf
u}_{g},{\bf x})p_{\theta}({\bf y}|{\bf u}_{l},{\bf u}_{g},{\bf z}),$ (1)

where $\theta$ denotes the collection of parameters for the generative model,
and $p_{\theta}({\bf u}_{l},{\bf u}_{g})=p_{\theta}({\bf
u}_{l})p_{\theta}({\bf u}_{g})$. Here $p_{\theta}({\bf z}|{\bf u}_{l},{\bf
u}_{g},{\bf x})$ and $p_{\theta}({\bf y}|{\bf u}_{l},{\bf u}_{g},{\bf z})$ are
the encoder and the predictor, respectively. Specifically, we have

$\displaystyle p_{\theta}({\bf u}_{l})$ $\displaystyle={\mathcal{N}}({\bf
0},{\bf I}),~{}~{}~{}~{}p_{\theta}({\bf u}_{g})={\mathcal{N}}({\bf 0},{\bf
I}),$ $\displaystyle p_{\theta}({\bf z}|{\bf u}_{l},{\bf u}_{g},{\bf x})$
$\displaystyle={\mathcal{N}}(\mu_{z}({\bf u}_{l},{\bf u}_{g},{\bf
x};\mbox{\boldmath$\theta$\unboldmath}),\sigma_{z}({\bf u}_{l},{\bf
u}_{g},{\bf x};\mbox{\boldmath$\theta$\unboldmath})),$ $\displaystyle
p_{\theta}({\bf y}|{\bf u}_{l},{\bf u}_{g},{\bf z})$
$\displaystyle={\mathcal{N}}(\mu_{y}({\bf u}_{l},{\bf u}_{g},{\bf
z};\mbox{\boldmath$\theta$\unboldmath}),\sigma_{y}({\bf u}_{l},{\bf
u}_{g},{\bf z};\mbox{\boldmath$\theta$\unboldmath})),$

where $\mu_{z}(\cdot;\cdot)$, $\sigma_{z}(\cdot;\cdot)$,
$\mu_{y}(\cdot;\cdot)$, $\sigma_{y}(\cdot;\cdot)$ are neural networks
parameterized by $\theta$. Note that for classification models the predictor
$p_{\theta}({\bf y}|{\bf u}_{l},{\bf u}_{g},{\bf z})$ becomes a categorical
distribution $Cat(f_{y}({\bf u}_{l},{\bf u}_{g},{\bf
z};\mbox{\boldmath$\theta$\unboldmath}))$, where $f_{y}(\cdot;\cdot)$ is a
neural network.

Inference Model for CounTS. We use an inference model $q_{\phi}({\bf y},{\bf
u}_{l},{\bf u}_{g},{\bf z}|{\bf x})$ to approximate the posterior distribution
of the latent variables, i.e., $p_{\theta}({\bf u}_{l},{\bf u}_{g},{\bf
z}|{\bf x})$. As shown in Fig. 1(right), we factorize $q_{\phi}({\bf y},{\bf
u}_{l},{\bf u}_{g},{\bf z}|{\bf x})$ as

$\displaystyle q_{\phi}({\bf y},{\bf u}_{l},{\bf u}_{g},{\bf z}|{\bf x})$
$\displaystyle=q_{\phi}({\bf y}|{\bf x})q_{\phi}({\bf u}_{l},{\bf u}_{g}|{\bf
x},{\bf y})q_{\phi}({\bf z}|{\bf x},{\bf y}),$ (2) $\displaystyle
q_{\phi}({\bf u}_{l},{\bf u}_{g}|{\bf x},{\bf y})$
$\displaystyle=q_{\phi}({\bf u}_{l}|{\bf x},{\bf y})q_{\phi}({\bf u}_{g}|{\bf
x},{\bf y})$ (3)

where $\phi$ is the collection of the inference model’s parameters, and . We
parameterized each factor in Eqn. 2 as

$\displaystyle q_{\phi}({\bf y}|{\bf x})$
$\displaystyle={\mathcal{N}}(\mu_{y}({\bf
x};\mbox{\boldmath$\phi$\unboldmath}),\sigma_{y}^{2}({\bf
x};\mbox{\boldmath$\phi$\unboldmath})),$ $\displaystyle q_{\phi}({\bf
u}_{l}|{\bf x},{\bf y})$ $\displaystyle={\mathcal{N}}(\mu_{u_{l}}({\bf x},{\bf
y};\mbox{\boldmath$\phi$\unboldmath}),\sigma_{u_{l}}^{2}({\bf x},{\bf
y};\mbox{\boldmath$\phi$\unboldmath})),$ $\displaystyle q_{\phi}({\bf
u}_{g}|{\bf x},{\bf y})$ $\displaystyle={\mathcal{N}}(\mu_{u_{g}}({\bf x},{\bf
y};\mbox{\boldmath$\phi$\unboldmath}),\sigma_{u_{g}}^{2}({\bf x},{\bf
y};\mbox{\boldmath$\phi$\unboldmath})),$ $\displaystyle q_{\phi}({\bf z}|{\bf
x},{\bf y})$ $\displaystyle={\mathcal{N}}(\mu_{z}({\bf x},{\bf
y};\mbox{\boldmath$\phi$\unboldmath}),\sigma_{z}^{2}({\bf x},{\bf
y};\mbox{\boldmath$\phi$\unboldmath})),$

where $\mu_{\cdot}(\cdot;\cdot)$ and $\sigma_{\cdot}(\cdot;\cdot)$ denote
neural networks with $\phi$ as their parameters.

Evidence Lower Bound. Our CounTS uses the evidence lower bound (ELBO)
${\mathcal{L}}_{ELBO}$ of the log likelihood $\log p({\bf y}|{\bf x})$ as an
objective to learn the generative and inference models. Maximizing the ELBO is
equivalent to learning the optimal variational distribution $q_{\phi}({\bf
u}_{l},{\bf u}_{g},{\bf z}|{\bf x})=\int q_{\phi}({\bf y},{\bf u}_{l},{\bf
u}_{g},{\bf z}|{\bf x})d{\bf y}$ that best approximates the posterior
distribution of the label and latent variables $p_{\theta}({\bf u}_{l},{\bf
u}_{g},{\bf z}|{\bf x})$. Specifically, we have

$\displaystyle{\mathcal{L}}_{ELBO}=$
$\displaystyle~{}{\mathbb{E}}_{q_{\phi}({\bf y},{\bf u}_{l},{\bf u}_{g},{\bf
z}|{\bf x})}[p_{\theta}({\bf y},{\bf u}_{l},{\bf u}_{g},{\bf z}|{\bf x})]$
$\displaystyle-{\mathbb{E}}_{q_{\phi}({\bf y},{\bf u}_{l},{\bf u}_{g},{\bf
z}|{\bf x})}[q_{\phi}({\bf y},{\bf u}_{l},{\bf u}_{g},{\bf z}|{\bf x})].$ (4)

Note that different from typical ELBOs, we explicitly involves ${\bf y}$ and
use $q_{\phi}({\bf y},{\bf u}_{l},{\bf u}_{g},{\bf z}|{\bf x})$ rather than
$q_{\phi}({\bf u}_{l},{\bf u}_{g},{\bf z}|{\bf x})$ (see Appendix A for
details); this is to expose the factor $q_{\phi}({\bf y}|{\bf x})$ in Eqn. 2
to allow for additional supervision on ${\bf y}$ (more details below). With
the factorization in Eqn. 1 and Eqn. 2, we can decompose Eqn. 4 as:

$\displaystyle{\mathcal{L}}$ ${}_{ELBO}=~{}{\mathbb{E}}_{q_{\phi}({\bf y}|{\bf
x})}{\mathbb{E}}_{q_{\phi}({\bf u}_{l},{\bf u}_{g},{\bf z}|{\bf x},{\bf
y})}[p_{\theta}({\bf y}|{\bf u}_{l},{\bf u}_{g},{\bf z})]$ (5)
$\displaystyle-{\mathbb{E}}_{q_{\phi}({\bf y},{\bf u}_{l},{\bf u}_{g}|{\bf
x})}\big{[}KL[q_{\phi}({\bf z}|{\bf x},{\bf y})||p_{\theta}({\bf z}|{\bf
u}_{l},{\bf u}_{g},{\bf x})]\big{]}$ (6)
$\displaystyle-{\mathbb{E}}_{q_{\phi}({\bf y}|{\bf x})}\big{[}KL[q_{\phi}({\bf
u}_{l},{\bf u}_{g}|{\bf x},{\bf y})||p_{\theta}({\bf u}_{l},{\bf
u}_{g})]\big{]}$ (7) $\displaystyle-{\mathbb{E}}_{q_{\phi}({\bf y}|{\bf
x})}[q_{\phi}({\bf y}|{\bf x})],$ (8)

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.

  1. (1)

Eqn. 5 predicts the label ${\bf y}$ using ${\bf u}_{l}$, ${\bf u}_{g}$, and
${\bf z}$ inferred from ${\bf x}$ (${\bf y}$ is marginalized out).

  2. (2)

Eqn. 6 regularizes the inference model $q_{\phi}({\bf z}|{\bf x},{\bf y})$ to
get closer to the generative model $p_{\theta}({\bf z}|{\bf u}_{l},{\bf
u}_{g},{\bf x})$.

  3. (3)

Eqn. 7 regularizes $q_{\phi}({\bf u}_{l},{\bf u}_{g}|{\bf x},{\bf y})$ using
the prior distribution $p({\bf u}_{l},{\bf u}_{g})$.

  4. (4)

Eqn. 8 regularizes $q_{\phi}({\bf y}|{\bf x})$ by maximizing its entropy,
preventing it from collapsing to deterministic solutions.

Figure 2: Network structure. We omit subscripts of $p_{\theta}$ and $q_{\phi}$
for clarity. Left: CounTS makes _predictions_ on ${\bf y}$ using the yellow
branch. Right: Given the current ${\bf x}$ and ${\bf y}^{pred}$, CounTS
generates _counterfactual explanation_ ${\bf x}^{cf}$ for any target label
${\bf y}^{cf}$ using the blue branch.

Final Objective Function. Inspired by (Louizos et al., 2017), we include an
additional term ($N$ is the training set size)

$\displaystyle{\mathcal{L}}_{y}=\sum\nolimits_{i=1}^{N}\log q({\bf y}_{i}|{\bf
x}_{i})$ (9)

to supervise $q({\bf y}|{\bf x})$ using label ${\bf y}$ during training; this
helps produce more accurate estimation for ${\bf u}_{l}$, ${\bf u}_{g}$, ${\bf
z}$ using $q_{\phi}({\bf u}_{l},{\bf u}_{g}|{\bf x},{\bf y})$ and
$q_{\phi}({\bf z}|{\bf x},{\bf y})$. The final objective function then becomes

$\displaystyle{\mathcal{L}}_{CounTS}={\mathcal{L}}_{ELBO}+\lambda{\mathcal{L}}_{y},$
(10)

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. 13, our self-
interpretable CounTS can perform inference to (1) make predictions on ${\bf
y}$ using the yellow branch in Fig. 2(left), and (2) generate counterfactual
explanation ${\bf x}^{cf}$ for any target label ${\bf 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_ ${\bf y}$. Specifically,
given an input time series ${\bf x}$, the encoder will first infer the
_initial label_ ${\bf y}$ using $q_{\phi}({\bf y}|{\bf x})$ and then infer
$({\bf u}_{l},{\bf u}_{g})$ and ${\bf z}$ from ${\bf x}$ and ${\bf y}$ using
$q_{\phi}({\bf u}_{l},{\bf u}_{g},{\bf z}|{\bf x},{\bf y})$. $({\bf
u}_{l},{\bf u}_{g})$ and ${\bf z}$ are then fed into the predictor
$p_{\theta}({\bf y}|{\bf u}_{l},{\bf u}_{g},{\bf z})$ to infer the _final_
label ${\bf y}$. Formally, CounTS predict the label as

$\displaystyle{\bf y}^{pred}={\mathbb{E}}_{q_{\phi}({\bf y}|{\bf
x})}{\mathbb{E}}_{q_{\phi}({\bf u}_{l},{\bf u}_{g},{\bf z}|{\bf x},{\bf
y})}[p_{\theta}({\bf y}|{\bf u}_{l},{\bf u}_{g},{\bf z})].$ (11)

Empirically, we find that directly using the means of $q_{\phi}({\bf y}|{\bf
x})$ and $q_{\phi}({\bf u}_{l},{\bf u}_{g},{\bf z}|{\bf x},{\bf y})$ as input
to $p_{\theta}({\bf y}|{\bf u}_{l},{\bf u}_{g},{\bf 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 ${\bf x}^{cf}$ defined below.

######  Definition 4.1 (Optimal Counterfactual Explanation).

Given a factual observation ${\bf x}$ and prediction ${\bf y}^{pred}$, the
optimal counterfactual explanation ${\bf x}^{cf}$ for the counterfactual
outcome for ${\bf y}^{cf}$ is

$\displaystyle{\bf x}^{cf}=\mathop{\rm argmax}\nolimits_{{\bf
x}^{\prime}}~{}p(Y_{{\bf x}={\bf x}^{\prime}}({\bf u})={\bf y}^{cf}),$

where ${\bf u}=({\bf u}_{l},{\bf u}_{g})$ and the counterfactual likelihood is
defined as

$\displaystyle p(Y_{{\bf x}={\bf x}^{\prime}}({\bf u})={\bf y}^{cf})$ (12)
$\displaystyle=\sum_{{\bf u}}p\left({\bf y}={\bf y}^{cf}|do\left({\bf x}={\bf
x}^{\prime}\right),{\bf u}\right)p({\bf u}|{\bf x}={\bf x},{\bf y}={\bf
y}^{pred}).$

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

######  Theorem 4.1 (Identifiability).

Given the posterior distribution of exogenous variable $p({\bf u}_{l},{\bf
u}_{g}|{\bf x},{\bf y})$, the effect of action $p({\bf y}={\bf y}^{cf}|do({\bf
x}={\bf x}^{\prime}),{\bf u}_{l},{\bf u}_{g})$ can be identified using
${\mathbb{E}}_{p({\bf z}|{\bf x}^{\prime},{\bf u}_{l},{\bf
u}_{g})}[p\left({\bf y}^{cf}|{\bf z},{\bf u}_{l},{\bf u}_{g}\right)]$.

See Appendix B.1 for the proof. With Theorem 4.1, we can rewrite Eqn. 12 as

$\displaystyle{\mathcal{L}}_{cf}={\mathbb{E}}_{p({\bf u}|{\bf x}={\bf x},{\bf
y}={\bf y}^{pred})}{\mathbb{E}}_{p({\bf z}|{\bf x}^{\prime},{\bf u})}[p({\bf
y}^{cf}|{\bf z},{\bf u})],$ (13)

where ${\bf u}=({\bf u}_{l},{\bf u}_{g})$ and $p({\bf u}|{\bf x}={\bf x},{\bf
y}={\bf y}^{pred})$ is approximated by $q_{\phi}({\bf u}|{\bf x}={\bf x},{\bf
y}={\bf y}^{pred})$. We use Monte Carlo estimates to compute the expectation
in Eqn. 12 and Eqn. 13, iteratively compute the gradient
$\frac{\partial{\mathcal{L}}_{cf}}{\partial{\bf x}^{\prime}}$ (via back-
propagation) to search for the optimal ${\bf x}^{\prime}$ in a way similar to
(Wang et al., 2019a; Mao et al., 2021a), and use it as ${\bf x}^{cf}$ (see the
complete algorithm in Appendix B.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) (Wachter et al., 2017a), Gradient-weighted Class Activation Mapping
(GradCAM) (Selvaraju et al., 2017), Gradient SHapley Additive exPlanations
(GradSHAP) (Lundberg & Lee, 2017a), Local Interpretable Model-agnostic
Explanations (LIME) (Ribeiro et al., 2016), Feature Importance in Time (FIT)
(Tonekaboni et al., 2020), Case-crossover APriori (CAP) (Dhaou et al., 2021),
and Counterfactual Residual Generative Adversarial Network (CounteRGAN)
(Nemirovsky et al., 2022) (see Appendix C.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 (Paszke et al., 2019). 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.

Table 1: Results on the toy dataset. We mark the best result with bold face
and the second best results with underline.

| CounteRGAN | RGD | GradCAM | GradSHAP | LIME | FIT | CAP | CounTS (Ours)
---|---|---|---|---|---|---|---|---
Pred. Accuracy (%) $\uparrow$ | ————————– 85.19 ————————– |  83.26
CCR $\uparrow$ |  1.25 | 1.21 | 1.09 | 1.13 | 1.2 | 1.15 | 0.97 | 1.33
Counterf. Accuracy (%) $\uparrow$ | 78.75 |  77.96 | - | - | - | 45.62 | - | 70.28

###  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 ${\bf 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 ${\bf x}$
affects the label, we introduce a mask vector ${\bf m}\in{\mathbb{R}}^{12}$
with first $6$ entries ${\bf m}_{1:6}$ set to $1$ and the last 6 entries ${\bf
m}_{7:12}$ set to $0$. We then generate ${\bf z}\in{\mathbb{R}}^{12}$ as ${\bf
z}={\bf u}\cdot({\bf m}\odot{\bf x})$ and the label $y\sim Bern(\sigma({\bf
z}^{\top}{\bf 1}+u))$ where $\sigma$ is the sigmoid function. Here ${\bf x}$’s
first $6$ entries ${\bf x}_{1:6}$ is label-related and the last $6$ entries
${\bf x}_{7:12}$ is label-agnostic.

Evaluation Metrics. We use three evaluation metrics:

  * •

Prediction Accuracy. This is the percentage of time series correctly predicted
($y^{pred}=y$) in the test set.

  * •

Counterfactual Accuracy. For a prediction model $f$, the generated
counterfactual explanation ${\bf x}^{cf}$, and the target label $y^{cf}$,
counterfactual accuracy is the percentage of time series where ${\bf x}^{cf}$
successfully change the model’s prediction to $y^{cf}$ (i.e., $f({\bf
x}^{cf})=y^{cf}$).

  * •

Counterfactual Change Ratio (CCR). This measures how well the counterfactual
explanation ${\bf x}^{cf}$ changes the label-related input ${\bf x}_{1:6}$
while keeping the label-agnostic input ${\bf x}_{7:12}$ unchanged. Formally we
use the average ratio of $\frac{\|{\bf x}_{1:6}^{cf}-{\bf
x}_{1:6}\|_{1}}{\|{\bf x}_{7:12}^{cf}-{\bf x}_{7:12}\|_{1}}$ across the test
set.

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 ${\bf 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 ${\bf 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({\bf u}|{\bf x},{\bf 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.

Table 2: Results on the _Spike_ dataset. We mark the best result with bold
face and the second best results with underline.

|  | CounteRGAN | RGD | GradCAM | GradSHAP | LIME | FIT | CAP | CounTS (Ours)
---|---|---|---|---|---|---|---|---|---
Pred. MSE $\downarrow$ |  | ————————– 0.128 ————————– | 0.117
CCR $\uparrow$ | 2 Active | 2.206 |  2.217 | 2.043 | 1.989 | 1.872 | 1.863 | 1.614 | 2.322
1 Active |  0.705 | 0.683 | 0.654 | 0.615 | 0.473 | 0.55 | 0.497 | 0.730
Counterf. MSE $\downarrow$ |  | 0.074 | 0.074 | - | - | - | 0.394 | - |  0.103

Figure 3: Qualitative results on the _Spike_ dataset. Column 1 shows the
original input ${\bf x}$ for CounTS, RGD, and CounteRGAN, respectively. Column
2 shows the predictions ${\bf y}^{pred}$ and counterfactual (target) labels
${\bf y}^{cf}$. The prediction for RGD and CounteRGAN are identical since they
use the same prediction model. Column 3 shows the counterfactual input ${\bf
x}^{cf}$. Column 4 shows the changes on input, i.e., ${\bf x}^{cf}-{\bf x}$.
FIT _cannot_ provide actionable explanations (see the importance score
generated by FIT in Appendix D).

###  5.3 Spike Dataset

Dataset Description. Inspired by (Tonekaboni et al., 2020), we construct the
_Spike_ synthetic dataset. In the dataset, each time series contains $3$
channel, i.e., ${\bf 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 ${\bf 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 ${\bf 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 ${\bf 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 C.1).

Evaluation Metrics. We use three evaluation metrics:

  * •

Prediction MSE. We use the mean square error (MSE)
$\frac{1}{N}\sum_{i}^{N}\|{\bf y}_{i}^{pred}-{\bf y}_{i}\|_{2}^{2}$ to measure
the prediction error in the test set with $N$ time series.

  * •

Counterfactual MSE. Similar to Sec. 5.2, for a prediction model $f$, the
generated counterfactual explanation ${\bf x}^{cf}$, and the target label
${\bf y}^{cf}$, counterfactual MSE is defined as
$\frac{1}{N}\sum_{i}^{N}\|f({\bf x}_{i}^{cf}))-{\bf y}_{i}^{cf}\|_{2}^{2}$; it
measures how successfully ${\bf x}^{cf}$ changes the model’s prediction to
$y^{cf}$.

  * •

CCR. For a given input ${\bf x}_{i}$, we set the target counterfactual label
${\bf y}_{i}^{cf}$ by shifting ${\bf y}^{pred}$ by $20$ timestamps to the
right. If at timestamp $t$, there is a spike in an active channel in ${\bf
x}_{i}$ triggering the output ${\bf y}_{i}^{pred}$ to switch from $0$ to $1$,
an ideal counterfactual explanation ${\bf x}_{i}^{cf}$ should (1) suppress all
spikes between $[t,t+20)$, (2) create a new spike at $t+20$ timestamp in all
active channels of the original input ${\bf x}_{i}$, and (3) keep all inactive
channels unchanged. Therefore the counterfactual change ratio can be defined
as (with $N$ time series): $CCR=\frac{1}{N}\sum_{i=1}^{N}\frac{\|{\bf
m}_{i}\odot({\bf x}_{i}-{\bf x}_{i}^{cf})\|_{1}}{\|({\bf 1}-{\bf
m}_{i})\odot({\bf x}_{i}-{\bf x}_{i}^{cf})\|_{1}},$ where ${\bf m}_{i}=[{\bf
u},{\bf u},\dots,{\bf u}]\in{\mathbb{R}}^{3\times T}$ repeats ${\bf u}$ in
each time step to mask inactive channels.

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 ${\bf u}$
that decides whether a spike in a specific channel affects the label ${\bf
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 ${\bf u}$ (i.e., the mask) and
therefore have more flexibility to modify the input ${\bf 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 ${\bf y}$ from $0$ to $1$. Both
CounTS’s and RGD’s predictions ${\bf y}^{pred}$ (blue in Column 1) are very
close to the ground truth.

Following the evaluation protocol above, we set our counterfactual (target)
label ${\bf y}^{cf}$ by shifting ${\bf 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 ${\bf y}^{pred}$ (see
Appendix D Fig. 5 for details).

Table 3: Average CCR for both intra-dataset and cross-dataset settings. The
average CCR is calculated over 6 counterfactual action settings
(A$\rightarrow$D, L$\rightarrow$D, R$\rightarrow$D, R$\rightarrow$A,
D$\rightarrow$A, and L$\rightarrow$A) for every model and dataset setting. The
_last row_ are the average CCR over all 6 dataset settings for every method.

| CounteRGAN | RGD | GradCAM | GradSHAP | LIME | FIT | CAP | CounTS (Ours)
---|---|---|---|---|---|---|---|---
$\rightarrow$SOF |  1.153 | 1.135 | 1.126 | 0.908 | 0.932 | 1.072 | 1.034 | 1.313
$\rightarrow$SHHS | 1.231 |  1.241 | 1.056 | 0.975 | 1.020 | 1.193 | 0.929 | 1.390
$\rightarrow$MESA |  1.285 | 1.188 | 1.031 | 1.025 | 1.007 | 1.142 | 1.031 | 1.377
SOF |  1.126 | 1.088 | 0.887 | 0.996 | 1.078 | 1.053 | 0.963 | 1.281
SHHS |  1.133 | 1.111 | 0.890 | 1.068 | 1.059 | 1.067 | 1.000 | 1.333
MESA | 1.163 |  1.166 | 0.883 | 1.042 | 0.987 | 1.089 | 0.962 | 1.329
Average |  1.182 | 1.155 | 0.979 | 1.002 | 1.014 | 1.103 | 0.987 | 1.337

###  5.4 Real-World Datasets

Dataset Description. We evaluate our model on three real-world medical
datasets: Sleep Heart Health Study (_SHHS_) (Quan et al., 1997), Multi-ethnic
Study of Atherosclerosis (_MESA_) (Zhang et al., 2018a), and Study of
Osteoporotic Fractures (_SOF_) (Cummings et al., 1990), 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).

Table 4: CCR for each counterfactual action setting in _SHHS_.

| A$\rightarrow$D | L$\rightarrow$D | R$\rightarrow$D | R$\rightarrow$A | D$\rightarrow$A | L$\rightarrow$A | Average
---|---|---|---|---|---|---|---
CounteRGAN | 1.119 | 1.144 | 1.055 | 1.298 | 1.125 | 1.059 |  1.133
RGD | 1.006 | 1.130 | 0.985 | 1.180 |  1.155 | 1.201 | 1.111
GradCAM | 0.882 | 0.790 | 0.859 | 0.954 | 0.882 | 0.977 | 0.890
GradSHAP | 1.100 | 1.105 | 0.903 | 1.171 | 1.100 | 1.027 | 1.068
LIME | 0.855 | 1.178 | 1.067 | 1.125 | 0.855 |  1.273 | 1.059
FIT | 0.948 | 1.012 | 1.215 |  1.405 | 0.948 | 0.877 | 1.067
CAP | 0.885 | 0.829 |  1.236 | 0.967 | 0.885 | 1.199 | 1.000
CounTS (Ours) | 1.244 |  1.159 | 1.349 | 1.515 | 1.324 | 1.407 | 1.333

Corresponding to our causal graph in Fig. 1(left), ${\bf x}$ is the original
breathing signal, ${\bf z}$ is signal patterns (representation) learned from
the encoder, ${\bf u}$ can be ‘gender’, ‘age’, and even ‘dataset index’
(‘dataset index’ will be used during cross-dataset prediction; more details
later), and ${\bf 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:

  * •

Prediction Metric and Counterfactual Metric are similar to those in the toy
classification dataset.

  * •

CCR. Since ground-truth confounders are not available in real-world datasets,
we propose to compute CCR on two consecutive 30-second segments with different
sleep stages. For example, given the two segments $[{\bf x}_{i,1},{\bf
x}_{i,2}]$ with the corresponding two predicted sleep stages $[{\bf
y}_{i,1}^{pred},{\bf y}_{i,2}^{pred}]=[A,D]$, we set the counterfactual labels
$[{\bf y}_{i,1}^{cf},{\bf y}_{i,2}^{cf}]=[A,A]$; we refer to this setting as
D$\rightarrow$A. An ideal counterfactual explanation $[{\bf x}_{i,1}^{cf},{\bf
x}_{i,2}^{cf}]$ should shift the model prediction from $[A,D]$ to $[A,A]$
while keeping the first segment unchanged (i.e., smaller $\|{\bf x}_{i,1}-{\bf
x}_{i,1}^{cf}\|_{1}$). CCR can therefore be computed as (with $N$ segment
pairs in the test set): $CCR=\frac{1}{N}\sum_{i=1}^{N}\frac{\|{\bf
x}_{i,2}-{\bf x}_{i,2}^{cf}\|_{1}}{\|{\bf x}_{i,1}-{\bf x}_{i,1}^{cf}\|_{1}}.$

In the experiments, we focus on cases where the target label ${\bf
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 ${\bf u}$, which has causal effect on
both the learned representation ${\bf z}$ and the sleep stage ${\bf 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).

  * •

Intra-Dataset Setting. Training and test sets are from same dataset, and we do
not involve the dataset index as part of the exogenous variable (confounder)
${\bf u}$.

  * •

Cross-Dataset Setting. We choose two of the datasets as the source datasets,
with the remaining one as the target dataset (e.g., _MESA_ +_SHHS_
$\rightarrow$_SOF_). We use all source datasets (e.g., _SHHS_ and _MESA_) and
$10\%$ of the target dataset (e.g., _SOF_) as the training set and use the
remaining $90\%$ of the target dataset as the test set. We treat the dataset
index as part of ${\bf u}$, which is predicted during both training and
inference.

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 D.

Table 5: Prediction accuracy for both intra-dataset and cross-dataset
settings. All baselines (e.g., RGD, CounteRGAN, and FIT) explain the same
prediction model and therefore share the same prediction accuracy.
‘$\rightarrow$_SHHS_ ’ means ‘ _MESA_ +_SOF_ $\rightarrow$_SHHS_ ’ and
similarly for ‘$\rightarrow$_MESA_ ’ and ‘$\rightarrow$_SOF_ ’.

Method |  $\rightarrow$_SOF_ |  $\rightarrow$_SHHS_ |  $\rightarrow$_MESA_ | _SOF_ | _SHHS_ | _MESA_
---|---|---|---|---|---|---
All Baselines | 0.702 | 0.747 | 0.768 | 0.729 | 0.801 | 0.813
CounTS (Ours) | 0.719 | 0.753 | 0.741 | 0.732 | 0.789 | 0.799

Table 6: Average counterfactual accuracy for both intra-dataset and cross-
dataset settings. The average is calculated over all 6 counterfactual action
settings.

|  $\rightarrow$SOF |  $\rightarrow$SHHS |  $\rightarrow$MESA | SOF | SHHS | MESA | Average
---|---|---|---|---|---|---|---
RGD |  0.913 | 0.887 | 0.907 | 0.862 | 0.907 | 0.892 | 0.895
CounteRGAN | 0.897 |  0.868 | 0.880 | 0.844 | 0.888 | 0.868 | 0.874
FIT | 0.331 | 0.327 | 0.346 | 0.248 | 0.322 | 0.294 | 0.311
CounTS (Ours) | 0.916 | 0.887 |  0.889 |  0.852 |  0.903 |  0.887 |  0.889

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 ${\bf u}$ and hence enjoys more flexibility when
modifying the input ${\bf 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 ${\bf u}_{l}$ and ${\bf 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 D.

Figure 4: Example on real-world dataset from _MESA_ intra-dataset experiment.
For this example, ${\bf y}^{pred}$ is [D,A] and target ${\bf y}^{cf}$ is [D,D]
since we assume ‘Deep Sleep’ and ‘Active’ have the most different patterns.
FIT _cannot_ provide actionable explanations (see the importance score
generated by FIT in Appendix D).

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:

  1. 1.

The First 30-Second Segment: In the original breathing signal, the first
30-second segment ($0\sim 30$ seconds) is periodic and contains 7 breathing
cycles (i.e., 7 periods). Our CounTS managed to keep the first 30-second
segment signal mostly unchanged, maintaining its periodicity and 7 breathing
cycles. In contrast, RGD changes the first segment to the point that its
periodicity is mostly lost. CounteRGAN can maintain the first segment’s
periodicity; however it increases the number breathing cycles from 7 to 8,
making it infeasible.

  2. 2.

The Second 30-Second Segment: In the original breathing signal, the second
30-second segment ($30\sim 60$ seconds) is irregular since the patient is in
the ‘Awake’ sleep stage. CounTS managed to make the ‘Awake’ (i.e., ‘Active’)
signal much more periodic, successfully steering the prediction to ‘Deep
Sleep’. In contrast, RGD and CounteRGAN failed to make the second 30-second
segment ($30\sim 60$ seconds) more periodic.

  3. 3.

Capturing and Maintaining Global Temporal Patterns: Meanwhile, we can observe
that CounTS’s explanation in the second 30-second segment ($30\sim 60$
seconds) has a much more similar breathing frequency (measured by the number
of breathing cycles per minute) with the first 30-second segment ($0\sim 30$
seconds), compared to RGD and CounteRGAN. This observation shows that CounTS
is capable of capturing global temporal pattern (the breathing frequency of
the individual) and keeping such exogenous variables unchanged in its
explanation thanks to the global exogenous variable $\mathbf{u}_{g}$ in our
model.

##  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 (Mi
et al., 2022), 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 (Zhang et al., 2018b; Quan et al., 1998), MESA
(Zhang et al., 2018b; Chen et al., 2014) and SOF (Zhang et al., 2018b; Spira
et al., 2008).

## References

  * Chen et al. (2018) Chen, J., Song, L., Wainwright, M., and Jordan, M.  Learning to explain: An information-theoretic perspective on model interpretation.  In _International Conference on Machine Learning_ , pp. 883–892. PMLR, 2018.
  * Chen et al. (2014) Chen, X., Wang, R., Zee, P., Lutsey, P., Javaheri, S., Alcántara, C., Jackson, C., Williams, M., and Redline, S.  Racial/ethnic differences in sleep disturbances: The multi-ethnic study of atherosclerosis (mesa).  _Sleep_ , 38, 11 2014.  doi: 10.5665/sleep.4732.
  * Cummings et al. (1990) Cummings, S. R., Black, D. M., Nevitt, M. C., Browner, W. S., Cauley, J. A., Genant, H. K., Mascioli, S. R., Scott, J. C., Seeley, D. G., Steiger, P., et al.  Appendicular bone density and age predict hip fracture in women.  _JAMA_ , 263(5):665–668, 1990.
  * Dhaou et al. (2021) Dhaou, A., Bertoncello, A., Gourvénec, S., Garnier, J., and Le Pennec, E.  Causal and interpretable rules for time series analysis.  In _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_, pp. 2764–2772, 2021.
  * Ding et al. (2022) Ding, H., Ma, Y., Deoras, A., Wang, Y., and Wang, H.  Zero-shot recommender systems.  In _ICLR Workshop on Deep Generative Models for Highly Structured Data_ , 2022.
  * Ghorbani et al. (2019) Ghorbani, A., Wexler, J., Zou, J. Y., and Kim, B.  Towards automatic concept-based explanations.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Goyal et al. (2019) Goyal, Y., Wu, Z., Ernst, J., Batra, D., Parikh, D., and Lee, S.  Counterfactual visual explanations.  In _International Conference on Machine Learning_ , pp. 2376–2384. PMLR, 2019.
  * Gupta et al. (2021) Gupta, S., Wang, H., Lipton, Z., and Wang, Y.  Correcting exposure bias for link recommendation.  In _International Conference on Machine Learning_ , pp. 3953–3963. PMLR, 2021.
  * Heo et al. (2018) Heo, J., Lee, H. B., Kim, S., Lee, J., Kim, K. J., Yang, E., and Hwang, S. J.  Uncertainty-aware attention for reliable interpretation and prediction.  In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018.
  * Huang et al. (2019) Huang, H., Wang, H., and Mak, B.  Recurrent poisson process unit for speech recognition.  In _AAAI_ , volume 33, pp. 6538–6545, 2019.
  * Jin et al. (2022) Jin, X., Park, Y., Maddix, D., Wang, H., and Wang, Y.  Domain adaptation for time series forecasting via attention sharing.  In _International Conference on Machine Learning_ , pp. 10280–10297. PMLR, 2022.
  * Kingma & Welling (2013) Kingma, D. P. and Welling, M.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Lim et al. (2021) Lim, B., Arık, S. Ö., Loeff, N., and Pfister, T.  Temporal fusion transformers for interpretable multi-horizon time series forecasting.  _International Journal of Forecasting_ , 37(4):1748–1764, 2021.
  * Lin et al. (2022) Lin, W., Lan, H., Wang, H., and Li, B.  Orphicx: A causality-inspired latent variable model for interpreting graph neural networks.  In _CVPR_ , 2022.
  * Louizos et al. (2017) Louizos, C., Shalit, U., Mooij, J. M., Sontag, D., Zemel, R., and Welling, M.  Causal effect inference with deep latent-variable models.  _Advances in neural information processing systems_ , 30, 2017.
  * Lundberg & Lee (2017a) Lundberg, S. M. and Lee, S.-I.  A unified approach to interpreting model predictions.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017a.
  * Lundberg & Lee (2017b) Lundberg, S. M. and Lee, S.-I.  A unified approach to interpreting model predictions.  _Advances in neural information processing systems_ , 30, 2017b.
  * Mao et al. (2021a) Mao, C., Chiquier, M., Wang, H., Yang, J., and Vondrick, C.  Adversarial attacks are reversible with natural supervision.  In _ICCV_ , 2021a.
  * Mao et al. (2021b) Mao, C., Gupta, A., Cha, A., Wang, H., Yang, J., and Vondrick, C.  Generative interventions for causal learning.  In _CVPR_ , 2021b.
  * Mi et al. (2022) Mi, L., Wang, H., Tian, Y., and Shavit, N.  Training-free uncertainty estimation for neural networks.  In _AAAI_ , 2022.
  * Natesan Ramamurthy et al. (2020) Natesan Ramamurthy, K., Vinzamuri, B., Zhang, Y., and Dhurandhar, A.  Model agnostic multilevel explanations.  _Advances in neural information processing systems_ , 33:5968–5979, 2020.
  * Nemirovsky et al. (2022) Nemirovsky, D., Thiebaut, N., Xu, Y., and Gupta, A.  Countergan: Generating counterfactuals for real-time recourse and interpretability using residual gans.  In Cussens, J. and Zhang, K. (eds.), _Proceedings of the Thirty-Eighth Conference on Uncertainty in Artificial Intelligence_ , volume 180 of _Proceedings of Machine Learning Research_ , pp. 1488–1497. PMLR, 01–05 Aug 2022.
  * Pan et al. (2021) Pan, Q., Hu, W., and Chen, N.  Two birds with one stone: Series saliency for accurate and interpretable multivariate time series forecasting.  In _IJCAI_ , pp. 2884–2891, 2021.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. Curran Associates, Inc., 2019.
  * Pawlowski et al. (2020) Pawlowski, N., Coelho de Castro, D., and Glocker, B.  Deep structural causal models for tractable counterfactual inference.  _Advances in Neural Information Processing Systems_ , 33:857–869, 2020.
  * Pearl (2009) Pearl, J.  _Causality_.  Cambridge University Press, 2 edition, 2009.
  * Plumb et al. (2018) Plumb, G., Molitor, D., and Talwalkar, A. S.  Model agnostic supervised local explanations.  _Advances in neural information processing systems_ , 31, 2018.
  * Plumb et al. (2020) Plumb, G., Al-Shedivat, M., Cabrera, Á. A., Perer, A., Xing, E., and Talwalkar, A.  Regularizing black-box models for improved interpretability.  _Advances in Neural Information Processing Systems_ , 33:10526–10536, 2020.
  * Quan et al. (1998) Quan, S., Howard, B., Iber, C., Kiley, J., Nieto, F., O’Connor, G., Rapoport, D., Redline, S., Robbins, J., Samet, J., and Wahl, P.  The sleep heart health study: Design, rationale, and methods.  _Sleep_ , 20:1077–85, 01 1998.  doi: 10.1093/sleep/20.12.1077.
  * Quan et al. (1997) Quan, S. F., Howard, B. V., Iber, C., Kiley, J. P., Nieto, F. J., O’connor, G. T., Rapoport, D. M., Redline, S., Robbins, J., Samet, J. M., et al.  The sleep heart health study: design, rationale, and methods.  _Sleep_ , 20(12):1077–1085, 1997.
  * Rajapaksha & Bergmeir (2022) Rajapaksha, D. and Bergmeir, C.  Limref: Local interpretable model agnostic rule-based explanations for forecasting, with an application to electricity smart meter data.  _arXiv preprint arXiv:2202.07766_ , 2022.
  * Ribeiro et al. (2016) Ribeiro, M. T., Singh, S., and Guestrin, C.  ”why should i trust you?” explaining the predictions of any classifier.  In _Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining_ , pp. 1135–1144, 2016.
  * Selvaraju et al. (2017) Selvaraju, R. R., Cogswell, M., Das, A., Vedantam, R., Parikh, D., and Batra, D.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In _Proceedings of the IEEE international conference on computer vision_ , pp. 618–626, 2017.
  * Shrikumar et al. (2017) Shrikumar, A., Greenside, P., and Kundaje, A.  Learning important features through propagating activation differences.  In _International conference on machine learning_ , pp. 3145–3153. PMLR, 2017.
  * Spira et al. (2008) Spira, A. P., Blackwell, T., Stone, K. L., Redline, S., Cauley, J. A., Ancoli-Israel, S., and Yaffe, K.  Sleep-disordered breathing and cognition in older women.  _Journal of the American Geriatrics Society_ , 56(1):45–50, 2008.
  * Sundararajan et al. (2017) Sundararajan, M., Taly, A., and Yan, Q.  Axiomatic attribution for deep networks.  In _International conference on machine learning_ , pp. 3319–3328. PMLR, 2017.
  * Tonekaboni et al. (2020) Tonekaboni, S., Joshi, S., Campbell, K., Duvenaud, D. K., and Goldenberg, A.  What went wrong and when? instance-wise feature importance for time-series black-box models.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.), _Advances in Neural Information Processing Systems_ , volume 33, pp. 799–809. Curran Associates, Inc., 2020.
  * Wachter et al. (2017a) Wachter, S., Mittelstadt, B., and Russell, C.  Counterfactual explanations without opening the black box: Automated decisions and the gdpr.  _Harv. JL & Tech._, 31:841, 2017a.
  * Wachter et al. (2017b) Wachter, S., Mittelstadt, B. D., and Russell, C.  Counterfactual explanations without opening the black box: Automated decisions and the gdpr.  _Cybersecurity_ , 2017b.
  * Wang (2017) Wang, H.  _Bayesian Deep Learning for Integrated Intelligence: Bridging the Gap between Perception and Inference_.  PhD thesis, Hong Kong University of Science and Technology, 2017.
  * Wang & Yeung (2016) Wang, H. and Yeung, D.-Y.  Towards bayesian deep learning: A framework and some existing methods.  _TDKE_ , 28(12):3395–3408, 2016.
  * Wang & Yeung (2020) Wang, H. and Yeung, D.-Y.  A survey on bayesian deep learning.  _CSUR_ , 53(5):1–37, 2020.
  * Wang et al. (2015) Wang, H., Wang, N., and Yeung, D.  Collaborative deep learning for recommender systems.  In _KDD_ , pp. 1235–1244, 2015.
  * Wang et al. (2019a) Wang, H., Mao, C., He, H., Zhao, M., Jaakkola, T. S., and Katabi, D.  Bidirectional inference networks: A class of deep bayesian networks for health profiling.  In _AAAI_ , volume 33, pp. 766–773, 2019a.
  * Wang et al. (2019b) Wang, S., Zhou, T., and Bilmes, J.  Bias also matters: Bias attribution for deep neural network explanation.  In _International Conference on Machine Learning_ , pp. 6659–6667. PMLR, 2019b.
  * Wang et al. (2020) Wang, Y., Menkovski, V., Wang, H., Du, X., and Pechenizkiy, M.  Causal discovery from incomplete data: A deep learning approach.  2020\.
  * Weinberger et al. (2020) Weinberger, E., Janizek, J., and Lee, S.-I.  Learning deep attribution priors based on prior knowledge.  _Advances in Neural Information Processing Systems_ , 33:14034–14045, 2020.
  * Yang et al. (2022) Yang, Y., Yuan, Y., Zhang, G., Wang, H., Chen, Y.-C., Liu, Y., Tarolli, C. G., Crepeau, D., Bukartyk, J., Junna, M. R., et al.  Artificial intelligence-enabled detection and assessment of parkinson’s disease using nocturnal breathing signals.  _Nature medicine_ , 28(10):2207–2215, 2022.
  * Zhang et al. (2018a) Zhang, G.-Q., Cui, L., Mueller, R., Tao, S., Kim, M., Rueschman, M., Mariani, S., Mobley, D., and Redline, S.  The national sleep research resource: towards a sleep data commons.  _JAMA_ , 25(10):1351–1358, 2018a.
  * Zhang et al. (2018b) Zhang, G.-Q., Cui, L., Mueller, R., Tao, S., Kim, M., Rueschman, M., Mariani, S., Mobley, D., and Redline, S.  The national sleep research resource: Towards a sleep data commons.  _Journal of the American Medical Informatics Association_ , pp. 572–572, 08 2018b.  doi: 10.1145/3233547.3233725.
  * Zhao et al. (2021) Zhao, M., Hoti, K., Wang, H., Raghu, A., and Katabi, D.  Assessment of medication self-administration using artificial intelligence.  _Nature medicine_ , 27(4):727–735, 2021.

Supplementary Material

##  Appendix A Evidence Lower Bound

We can derive the evidence lower bound (ELBO) in Eqn. 4 as follows:

$\displaystyle\log p({\bf y}|{\bf x})$ $\displaystyle=\log\int_{{\bf
z}}\int_{{\bf u}}p({\bf y},{\bf z},{\bf u}|{\bf x})d{\bf u}d{\bf z}$
$\displaystyle\geq\mathbb{E}_{q_{\phi}({\bf u},{\bf z},{\bf y}|{\bf
x})}\left[\log\frac{p_{\theta}({\bf y},{\bf u},{\bf z}|{\bf x})}{q_{\phi}({\bf
u},{\bf z},{\bf y}|{\bf x})}\right]$ $\displaystyle=\mathbb{E}_{q_{\phi}({\bf
u},{\bf z},{\bf y}|{\bf x})}\left[\log\frac{p({\bf u})p_{\theta}({\bf z}|{\bf
u},{\bf x})p_{\theta}({\bf y}|{\bf u},{\bf z})}{q_{\phi}({\bf u},{\bf z},{\bf
y}|{\bf x})}\right]$ $\displaystyle={\mathbb{E}}_{q_{\phi}({\bf y},{\bf
u}_{l},{\bf u}_{g},{\bf z}|{\bf x})}[p_{\theta}({\bf y},{\bf u}_{l},{\bf
u}_{g},{\bf z}|{\bf x})]-{\mathbb{E}}_{q_{\phi}({\bf y},{\bf u}_{l},{\bf
u}_{g},{\bf z}|{\bf x})}[q_{\phi}({\bf y},{\bf u}_{l},{\bf u}_{g},{\bf z}|{\bf
x})]$ $\displaystyle=\mathbb{E}_{q_{\phi}({\bf u},{\bf z},{\bf y}|{\bf
x})}[\log p({\bf u})]+\mathbb{E}_{q_{\phi}({\bf u},{\bf z},{\bf y}|{\bf
x})}[\log p_{\theta}({\bf z}|{\bf u},{\bf x})]+\mathbb{E}_{q_{\phi}({\bf
u},{\bf z},{\bf y}|{\bf x})}[\log p_{\theta}({\bf y}|{\bf u},{\bf
z})]-\mathbb{E}_{q_{\phi}({\bf u},{\bf z},{\bf y}|{\bf x})}[\log q_{\phi}({\bf
u},{\bf z},{\bf y}|{\bf x})]$

##  Appendix B More Details on Counterfactual Inference

###  B.1 Proof on Identifiability

######  Theorem B.1 (Identifiability).

Given the posterior distribution of exogenous variable $p({\bf u}_{l},{\bf
u}_{g}|{\bf x},{\bf y})$, the effect of action $p({\bf y}={\bf y}^{cf}|do({\bf
x}={\bf x}^{\prime}),{\bf u}_{l},{\bf u}_{g})$ can be identified using
${\mathbb{E}}_{p({\bf z}|{\bf x}^{\prime},{\bf u}_{l},{\bf u}_{g})}p\left({\bf
y}^{cf}|{\bf z},{\bf u}_{l},{\bf u}_{g}\right)$.

###### Proof.

With ${\bf u}=({\bf u}_{l},{\bf u}_{g})$ and applying Rule 2 and 3 in do-
calculus (Pearl, 2009), we have

$\displaystyle p\left({\bf y}^{cf}|do\left({\bf x}={\bf
x}^{\prime}\right),{\bf u}\right)$ $\displaystyle\ \ =\int_{{\bf
z}}p\left({\bf y}^{cf}|{\bf z},do\left({\bf x}^{\prime}\right),{\bf
u}\right)p\left({\bf z}|do\left({\bf x}^{\prime}\right),{\bf u}\right)d{\bf
z}$ $\displaystyle\ \ ={\mathbb{E}}_{p\left({\bf z}|do\left({\bf
x}^{\prime}\right),{\bf u}\right)}\left[p\left({\bf y}^{cf}|{\bf
z},do\left({\bf x}^{\prime}\right),{\bf u}\right)\right]$
$\displaystyle\overset{\mathrm{Rule2}}{=}{\mathbb{E}}_{p\left({\bf z}|{\bf
x}^{\prime},{\bf u}\right)}\left[p\left({\bf y}^{cf}|do({\bf z}),do\left({\bf
x}^{\prime}\right),{\bf u}\right)\right]$
$\displaystyle\overset{\mathrm{Rule3}}{=}{\mathbb{E}}_{p\left({\bf z}|{\bf
x}^{\prime},{\bf u}\right)}\left[p\left({\bf y}^{cf}|do({\bf z}),{\bf
u}\right)\right]$
$\displaystyle\overset{\mathrm{Rule2}}{=}{\mathbb{E}}_{p\left({\bf z}|{\bf
x}^{\prime},{\bf u}\right)}\left[p\left({\bf y}^{cf}|{\bf z},{\bf
u}\right)\right],$ (14)

concluding the proof. ∎

###  B.2 Counterfactual Explanation Algorithm

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

Algorithm 1 Generating Counterfactual Explanations

Input: data ${\bf x}$, threshold $\epsilon$, number of samples $n$ and $m$.

Sample ${\bf y}^{pred}$ from $q_{\phi}({\bf y}|{\bf x})$.

Set a target counterfactual output ${\bf y}^{cf}\neq{\bf y}^{pred}$

for $i=1$ to $m$ do

Sample ${\bf u}_{i}$ from $q_{\phi}({\bf u}|{\bf x},{\bf y})$

for $j=1$ to $n$ do

Sample ${\bf z}_{ij}\sim p_{\theta}({\bf z}|{\bf x},{\bf u}_{i})$

Sample ${\bf y}_{ij}\sim p_{\theta}({\bf y}|{\bf z}_{ij},{\bf u}_{i})$

end for

end for

Calculate ${\bf y}^{cf}=\frac{\sum_{i}\sum_{j}{\bf y}_{ij}}{n\times m}$

Set ${\bf x}^{cf}={\bf x}$

while $|{\bf y}-{\bf y}^{cf}|\geq\epsilon$ do

Update ${\bf x}^{cf}={\bf x}^{cf}-\lambda\frac{\partial\mid{\bf y}-{\bf
y}^{cf}\mid}{\partial{\bf x}^{cf}}$

end while

return ${\bf x}^{cf}$

##  Appendix C More Details on Experiments

###  C.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$
(15)

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.

$\displaystyle\mathbf{\mbox{\boldmath$\theta$\unboldmath}}=[0.8,0.4,0];$ (16)
$\displaystyle n_{d}\sim\operatorname{Bernoulli}(\mathbf{\theta}_{d});$
$\displaystyle m_{d}\sim\operatorname{Bernoulli}(\mathbf{\theta}_{d});$
$\displaystyle\eta_{d}=\begin{cases}\operatorname{Poisson}(\lambda=2)&\mathbf{1}\text{
if }\left(n_{d}==1\right)\\\ 0&\text{ otherwise }\end{cases}$
$\displaystyle\mathbf{g}_{d}=\operatorname{Sample}\left([T],\eta_{d}\right)$
$\displaystyle 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}$
$\displaystyle
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.$

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).

###  C.2 Details on Baselines

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

  * •

Gradient-based methods. Regularized Gradient Descent (RGD) (Wachter et al.,
2017a) 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) (Selvaraju et al., 2017) 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) (Lundberg & Lee, 2017a) 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) (Ribeiro et al., 2016) 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) (Tonekaboni et
al., 2020) 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) (Dhaou et al., 2021)
applies association rule mining algorithm, Apriori, to explore the causal
relationship in time-series data.

  * •

Generative method. Counterfactual Residual Generative Adversarial Network
(CounteRGAN) (Nemirovsky et al., 2022) 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.

##  Appendix D Additional Results

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

Figure 5: Qualitative results from FIT in the _Spike_ dataset. Figure 6:
Qualitative results on FIT in a _MESA_ intra-dataset setting.

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. 6.

Table 7$\sim$11 show more detailed CCR results with different counterfactual
action settings and dataset settings. Table 12 shows detailed counterfactual
accuracy for each dataset settings and each counterfactual action settings.

Table 7: CCR for different counterfactual action settings (e.g.,
A$\rightarrow$D means ‘Awake’$\rightarrow$‘Deep Sleep’) in _SOF_.

| A$\rightarrow$D | L$\rightarrow$D | R$\rightarrow$D | R$\rightarrow$A | D$\rightarrow$A | L$\rightarrow$A | Average
---|---|---|---|---|---|---|---
CounteRGAN | 1.356 |  1.138 | 1.165 | 1.127 | 0.908 | 1.060 |  1.126
RGD | 0.996 | 1.103 |  1.174 | 0.951 | 1.108 |  1.194 | 1.088
GradCAM | 1.012 | 0.926 | 0.644 | 0.727 | 1.012 | 1.003 | 0.887
GradSHAP | 0.843 | 1.177 |  1.219 | 0.913 | 0.843 | 0.981 | 0.996
LIME |  1.220 | 0.857 | 1.095 | 0.936 |  1.220 | 1.139 | 1.078
FIT | 1.064 | 1.150 | 1.169 |  1.202 | 1.064 | 0.667 | 1.053
CAP | 0.989 | 0.952 | 1.019 | 1.139 | 0.989 | 0.689 | 0.963
CounTS (Ours) | 1.179 | 1.161 | 1.240 | 1.475 | 1.245 | 1.387 | 1.281

Table 8: CCR for different counterfactual action settings in _SHHS+SOF
$\rightarrow$ MESA_.

| A$\rightarrow$D | L$\rightarrow$D | R$\rightarrow$D | R$\rightarrow$A | D$\rightarrow$A | L$\rightarrow$A | Average
---|---|---|---|---|---|---|---
CounteRGAN |  1.381 | 1.379 | 1.068 | 1.373 |  1.355 | 1.158 |  1.285
RGD | 1.226 | 1.130 | 1.089 | 1.215 | 1.177 |  1.293 | 1.188
GradCAM | 1.056 | 1.014 | 0.743 | 1.218 | 1.056 | 1.101 | 1.031
GradSHAP | 1.174 | 0.822 | 1.046 | 0.948 | 1.174 | 0.988 | 1.025
LIME | 0.904 | 1.246 | 1.138 | 1.017 | 0.904 | 0.836 | 1.007
FIT | 1.143 | 1.228 |  1.193 | 1.206 | 1.143 | 0.938 | 1.142
CAP | 0.969 | 1.179 | 0.947 |  1.418 | 0.969 | 0.703 | 1.031
CounTS (Ours) | 1.391 |  1.331 | 1.253 | 1.470 | 1.465 | 1.355 | 1.377

Table 9: CCR for different counterfactual action settings in _MESA_.

| A$\rightarrow$D | L$\rightarrow$D | R$\rightarrow$D | R$\rightarrow$A | D$\rightarrow$A | L$\rightarrow$A | Average
---|---|---|---|---|---|---|---
CounteRGAN |  1.183 | 1.345 | 1.160 | 1.117 | 0.974 | 1.200 | 1.163
RGD | 1.159 | 1.088 | 1.149 | 1.204 | 1.190 |  1.207 |  1.166
GradCAM | 1.078 | 0.790 | 0.419 | 0.954 | 1.078 | 0.977 | 0.883
GradSHAP | 0.885 | 1.027 | 1.039 |  1.273 | 0.885 | 1.145 | 1.042
LIME | 0.903 | 1.003 | 0.930 | 1.233 | 0.903 | 0.952 | 0.987
FIT | 1.210 | 1.118 | 1.024 | 1.082 |  1.210 | 0.887 | 1.089
CAP | 0.786 | 1.001 |  1.176 | 1.243 | 0.786 | 0.783 | 0.962
CounTS (Ours) | 1.213 |  1.309 | 1.294 | 1.526 | 1.350 | 1.285 | 1.329

Table 10: CCR for different counterfactual action settings in _MESA+SOF
$\rightarrow$SHHS_.

| A$\rightarrow$D | L$\rightarrow$D | R$\rightarrow$D | R$\rightarrow$A | D$\rightarrow$A | L$\rightarrow$A | Average
---|---|---|---|---|---|---|---
CounteRGAN | 1.317 | 1.157 | 1.281 | 1.300 | 1.270 | 1.060 | 1.231
RGD | 1.193 |  1.266 | 1.115 |  1.313 |  1.304 |  1.255 |  1.241
GradCAM | 1.131 | 1.093 | 0.725 | 1.306 | 1.131 | 0.953 | 1.056
GradSHAP | 1.025 | 0.910 | 1.092 | 0.886 | 1.025 | 0.910 | 0.975
LIME | 1.018 | 1.092 | 1.118 | 0.944 | 1.018 | 0.930 | 1.020
FIT | 1.270 | 1.259 | 1.105 | 1.263 | 1.270 | 0.991 | 1.193
CAP | 0.994 | 0.819 | 1.033 | 0.953 | 0.994 | 0.777 | 0.929
CounTS (Ours) |  1.301 | 1.452 |  1.167 | 1.560 | 1.417 | 1.443 | 1.390

Table 11: CCR for different counterfactual action settings in _SHHS+MESA
$\rightarrow$SOF_.

| A$\rightarrow$D | L$\rightarrow$D | R$\rightarrow$D | R$\rightarrow$A | D$\rightarrow$A | L$\rightarrow$A | Average
---|---|---|---|---|---|---|---
CounteRGAN | 1.210 | 1.050 | 1.034 | 1.403 | 1.065 | 1.248 |  1.153
RGD | 1.183 | 0.993 | 1.015 | 1.260 |  1.155 | 1.201 | 1.135
GradCAM | 1.108 | 1.018 | 1.020 | 1.215 | 1.108 |  1.286 | 1.126
GradSHAP | 0.983 | 0.929 | 0.780 | 1.044 | 0.983 | 0.728 | 0.908
LIME | 0.738 | 1.051 |  1.199 | 1.155 | 0.738 | 0.708 | 0.932
FIT | 1.110 |  1.175 | 0.954 | 1.206 | 1.110 | 0.876 | 1.072
CAP | 1.106 | 0.951 | 0.833 | 1.145 | 1.106 | 1.060 | 1.034
CounTS (Ours) |  1.203 | 1.266 | 1.202 |  1.395 | 1.472 | 1.341 | 1.313

Table 12: Couterfactual accuracy for each dataset settings and each
counterfactual action settings.

| Method | A$\rightarrow$D | L$\rightarrow$D | R$\rightarrow$D | R$\rightarrow$A | D$\rightarrow$A | L$\rightarrow$A | Average
---|---|---|---|---|---|---|---|---
SHHS+MESA $\downarrow$ SOF | RGD |  0.867 |  0.903 |  0.921 | 0.932 | 0.931 | 0.922 |  0.913
CounteRGAN | 0.848 | 0.917 | 0.903 | 0.917 | 0.896 | 0.900 | 0.897
FIT | 0.285 | 0.296 | 0.328 | 0.384 | 0.285 | 0.348 | 0.331
CounTS (Ours) | 0.868 | 0.893 | 0.929 |  0.918 |  0.920 |  0.913 | 0.916
MESA+SOF $\downarrow$ SHHS | RGD |  0.828 | 0.882 | 0.916 | 0.910 | 0.903 | 0.890 | 0.887
CounteRGAN |  0.828 | 0.869 |  0.920 | 0.867 | 0.880 |  0.875 |  0.868
FIT | 0.334 | 0.263 | 0.374 | 0.354 | 0.334 | 0.348 | 0.327
CounTS (Ours) | 0.837 |  0.878 | 0.923 |  0.899 |  0.899 | 0.890 | 0.887
SHHS+SOF $\downarrow$ MESA | RGD | 0.882 | 0.898 |  0.933 | 0.910 | 0.938 |  0.888 | 0.907
CounteRGAN |  0.872 | 0.878 | 0.926 | 0.905 | 0.901 | 0.897 | 0.880
FIT | 0.335 | 0.269 | 0.301 | 0.357 | 0.335 | 0.397 | 0.346
CounTS (Ours) | 0.882 |  0.888 | 0.936 | 0.910 |  0.930 | 0.881 |  0.889
SOF | RGD |  0.823 |  0.820 | 0.919 |  0.856 |  0.861 | 0.891 | 0.862
CounteRGAN | 0.831 | 0.800 | 0.899 | 0.867 | 0.832 | 0.875 | 0.844
FIT | 0.257 | 0.328 | 0.311 | 0.328 | 0.257 | 0.310 | 0.248
CounTS (Ours) | 0.813 | 0.831 | 0.907 | 0.843 | 0.864 |  0.881 |  0.852
SHHS | RGD | 0.880 | 0.892 | 0.940 | 0.944 |  0.887 | 0.900 | 0.907
CounteRGAN | 0.841 | 0.854 | 0.916 | 0.913 | 0.888 | 0.875 | 0.888
FIT | 0.335 | 0.272 | 0.353 | 0.337 | 0.335 | 0.348 | 0.322
CounTS (Ours) |  0.868 |  0.880 |  0.930 | 0.931 | 0.875 |  0.889 |  0.903
MESA | RGD | 0.846 | 0.888 |  0.912 | 0.907 | 0.909 |  0.888 | 0.892
CounteRGAN | 0.828 | 0.869 | 0.920 | 0.867 | 0.880 | 0.875 | 0.868
FIT | 0.290 | 0.316 | 0.320 | 0.307 | 0.290 | 0.301 | 0.294
CounTS (Ours) |  0.837 |  0.878 | 0.903 |  0.899 |  0.899 | 0.890 |  0.887