Counterfactual Explanations in Sequential Decision Making Under Uncertainty

Methods to find counterfactual explanations have predominantly focused on one step decision making processes. In this work, we initiate the development of methods to find counterfactual explanations for decision making processes in which multiple, dependent actions are taken sequentially over time. We start by formally characterizing a sequence of actions and states using finite horizon Markov decision processes and the Gumbel-Max structural causal model. Building upon this characterization, we formally state the problem of finding counterfactual explanations for sequential decision making processes. In our problem formulation, the counterfactual explanation specifies an alternative sequence of actions differing in at most k actions from the observed sequence that could have led the observed process realization to a better outcome. Then, we introduce a polynomial time algorithm based on dynamic programming to build a counterfactual policy that is guaranteed to always provide the optimal counterfactual explanation on every possible realization of the counterfactual environment dynamics. We validate our algorithm using both synthetic and real data from cognitive behavioral therapy and show that the counterfactual explanations our algorithm finds can provide valuable insights to enhance sequential decision making under uncertainty.


Introduction
In recent years, there has been a rising interest on the potential of machine learning models to assist and enhance decision making in high-stakes applications such as justice, education and healthcare [1][2][3]. In this context, machine learning models and algorithms have been used not only to predict the outcome of a decision making process from a set of observable features but also to find what would have to change in (some of) the features for the outcome of a specific process realization to change. For example, in loan decisions, a bank may use machine learning both to estimate the probability that a customer repays a loan and to find by how much a customer may need to reduce her credit card debt to increase the probability of repayment over a given threshold. Here, our focus is on the latter, which has been often referred to as counterfactual explanations.
The rapidly growing literature on counterfactual explanations has predominantly focused on one-step decision making processes as the one described above [4][5][6]. In such settings, the probability that an outcome occurs is typically estimated using a supervised learning model and finding counterfactual explanations reduces to a search problem across the space of features and model predictions [7][8][9][10][11]. Moreover, it has been argued that, to obtain counterfactual explanations that are actionable and have the predicted effect on the outcome, one should favor causal models [12][13][14]. In this work, our goal is instead to find counterfactual explanations for decision making processes in which multiple, dependent actions are taken sequentially over time. In this setting, the (final) outcome of the process depends on the overall sequence of actions and the counterfactual explanation specifies an alternative sequence of actions differing in at most k actions from the observed sequence that could have led the process realization to a better outcome. For example, in medical treatment, assume a physician takes a sequence of actions to treat a patient but the patient's prognosis does not improve, then a counterfactual explanation would help the physician understand how a small number of actions taken differently could have improved the patient's prognosis. However, since there is (typically) uncertainty on the counterfactual dynamics of the environment, there may be a different counterfactual explanation that is optimal under each possible realization of the counterfactual dynamics. Moreover, since, in many realistic scenarios, a decision maker needs to decide among a small number of actions on the basis of a few observable covariates, which are often discretized into (percentile) ranges, in this work, we focus on decision making processes where the state and action spaces are discrete and low-dimensional.
The work most closely related to ours, which lies within the area of machine learning for healthcare, has achieved early success at using machine learning to enhance sequential (clinical) decisions [15][16][17][18]. However, rather than finding counterfactual explanations, this line of work has predominantly focused on using reinforcement learning to design better treatment policies. A notable exception is by Oberst and Sontag [18], which introduces an off-policy evaluation procedure for highlighting specific realizations of a sequential decision making process where a policy trained using reinforcement learning would have led to a substantially different outcome. While our work builds upon their modeling framework, we do not focus on off-policy evaluation but, instead, on the generation of (better) alternative action sequences as a means of explanation for the process's outcome. Therefore, we see their contributions as complementary to ours.
More broadly, our work is not the first to use counterfactual reasoning in the context of sequential decision making in the machine learning literature [19,20]. However, previous work has used counterfactual reasoning either to learn optimal policies using limited observational data [19] or to explain the action choices of a reinforcement learning agent [20]. In contrast, our work uses a causal model of the environment to compute alternative action sequences, close to the observed one, which, in retrospect, could have improved the outcome of the decision making process. Refer to Appendix A for a discussion of further related work.
Our contributions. We start by formally characterizing a sequence of (discrete) actions and (discrete) states using finite horizon Markov decision processes (MDPs). Here, we model the transition probabilities between a pair of states, given an action, using the Gumbel-Max structural causal model [18]. This model has been shown to have a desirable counterfactual stability property and, given a sequence of actions and states, it allows us to reliably estimate the counterfactual outcome under an alternative sequence of actions. Building upon this characterization, we formally state the problem of finding the counterfactual explanation for an observed realization of a sequential decision making process as a constrained search problem over the set of alternative sequences of actions differing in at most k actions from the observed sequence. Then, we present a polynomial time algorithm based on dynamic programming that finds a counterfactual policy that is guaranteed to always provide the optimal counterfactual explanation on every possible realization of the counterfactual transition probability induced by the observed process realization. Finally, we validate our algorithm using both synthetic and real data from cognitive behavioral therapy and show that the counterfactual explanations can provide valuable insights to enhance sequential decision making under uncertainty 1 .

Characterizing Sequential Decision Making using Causal Markov Decision Processes
Our starting point is the following stylized setting, which resembles a variety of real-world sequential decision making processes. At each time step t ∈ {0, . . . , T − 1}, the decision making process is characterized by a state s t ∈ S, where S is a space of n states, an action a t ∈ A, where A is a space of m actions, and a reward R(s t , a t ) ∈ R. Moreover, given a realization of a decision making process τ = {(s t , a t )} T −1 t=0 , the outcome of the decision making process o(τ ) = t R(s t , a t ) is given by the sum of the rewards.
Given the above setting, we characterize the relationship between actions, states and outcomes using finite horizon Markov decision processes (MDPs). More specifically, we consider an MDP M = (S, A, P, R, T ), where S is the state space, A is the set of actions, P denotes the transition probability P (S t+1 = s t+1 | S t = s t , A t = a t ), R denotes the immediate reward R(s t , a t ), and T is the time horizon. While this characterization Figure 1: Structural causal model C for a Markov decision process M. Green boxes represent endogenous random variables and pink boxes represent exogenous noise variables. The value of each endogenous variable is given by a function of the values of its ancestors in the structural causal model, as defined by Eq. 1. The value of each exogenous noise variable is sampled independently from a given distribution. An intervention do(A t = a) breaks the dependence of the variable A t from its ancestors (highlighted by dotted lines) and sets its value to a. After observing an event S t+1 = s t+1 , S t = s t , A t = a t , a counterfactual prediction can be thought of as the result of an intervention do(A t = a) in a modified SCM where U t takes values u t from a posterior distribution with support such that s t+1 = g S (s t , a t , u t ).
is helpful to make predictions about future states and design action policies [21], it is not sufficient to make counterfactual predictions, e.g., given a realization of a decision making process τ = {(s t , a t )} T −1 t=0 , we cannot know what would have happened if, instead of taking action a t at time t, we had taken action a = a t . To be able to overcome this limitation, we will now augment the above characterization using a particular class of structural causal model (SCM) [22,23] with desirable properties.
Let C be a structural causal model defined by the assignments where U t and V t are n-and m-dimensional independent noise variables, respectively, and g S and g A are two given functions, and let P C be the distribution entailed by C. Then, as argued by Buesing et al. [19], we can always find a distribution for the noise variables and a function g S so that the transition probability of the MDP of interest is given by the following interventional distribution over the SCM C: where do(A t = a t ) denotes a (hard) intervention in which the second assignment in Eq. 1 is replaced by the value a t . Under this view, given an observed realization of a decision making process τ = {(s t , a t )} T −1 t=0 , we can compute the posterior distribution P C | St=st,St+1=st+1,At=at (U t ) of each noise variable U t . Building on the conditional density function of this posterior distribution, which we denote as f C | St=st,St+1=st+1,At=at Ut (u), we can define a (non-stationary) counterfactual transition probability where, in (a), we drop the do(·) because U t and A t are independent in the modified SCM. Importantly, the above counterfactual transition probability allows us to make counterfactual predictions, e.g., given that, at time t the state was s t and, at time t + 1, the process transitioned to state s t+1 after taking action a t , what would have been the probability of transitioning to state s after taking action a = a t if the state had been s at time t. Refer to Figure 1 for a visual description of our model and the notion of counterfactual predictions. However, for state variables taking discrete values, the posterior distribution of the noise may be nonidentifiable without further assumptions, as argued by Oberst and Sontag [18]. This is because there may be several noise distributions and functions g S which give interventional distributions consistent with the MDP's transition probabilities but result in different counterfactual transition distributions. To avoid these non-identifiability issues, we follow Oberst and Sontag and restrict our attention to the class of Gumbel-Max SCMs, i.e., where U t,s ∼ Gumbel(0, 1) and P (· | S t , A t ) is the transition probability of the MDP. More specifically, this class of SCMs has been shown to satisfy a desirable counterfactual stability property, which goes intuitively as follows. Assume that, at time t, the process transitioned from state s t to state s t+1 after taking action a t . Then, in a counterfactual scenario, it is unlikely that, at time t, the process would transition from a state s to a state s other than s t+1 -the factual one-unless choosing an action a that decreases the relative chances of s t+1 compared to the other states. More formally, for a model satisfying counterfactual stability, given it holds that P τ,t (S t+1 = s | S t = s, A t = a) = 0. In practice, in addition to solving the non-identifiability issues, the use of Gumbel-Max SCMs allows for an efficient procedure to sample from the corresponding noise posterior distribution P C | St=st,St+1=st+1,At=at (U t ), described elsewhere [18,24], and given a set of d samples from the noise posterior distribution, we can compute an unbiased finite sample Monte-Carlo estimator for the counterfactual transition probability, as defined in Eq. 3, as follows:

Counterfactual Explanations in Markov Decision Processes
Inspired by previous work on counterfactual explanations in supervised learning [7,8], we focus on counterfactual outcomes that could have occurred if the alternative action sequence was "close" to the observed one. However, since in our setting, there is uncertainty on the counterfactual dynamics of the environment, we will look for a non-stationary counterfactual policy π that, under every possible realization of the counterfactual transition probability defined in Eq. 3, is guaranteed to provide the optimal alternative sequence of actions differing in at most k actions from the observed one. More specifically, let τ = {(s t , a t )} T −1 t=0 be an observed realization of a decision making process characterized by a Markov decision process (MDP) M = (S, A, P, R, T ) with a transition probability defined via a Gumbel-Max structural causal model (SCM), as described in Section 2. Then, to characterize the effect that any alternative action sequence would have had on the outcome of the above process realization, we start by building a non-stationary counterfactual MDP M τ = (S + , A, P + τ , R + , T ). Here, S + = S ×{0, . . . , T −1} is an enhanced state space such that each s + ∈ S + corresponds to a pair (s, l) indicating that the original decision making process would have been at state s ∈ S had already taken l actions differently from the observed sequence. Following this definition, R + denotes the reward function which we define as r + ((s, l), a) = R(s, a) for any (s, l) ∈ S + and a ∈ A, i.e., the counterfactual rewards remain independent of the number of modifications in the action sequence. Lastly, let P τ be the counterfactual transition probability, as defined by Eq. 3. Then, the transition probability P + τ for the enhanced state space is defined as: if (a = a t ∧ l = l) ∨ (a = a t ∧ l = l + 1) 0 otherwise, ALGORITHM 1: It samples a counterfactual explanation from the counterfactual policy π Input: counterfactual policy π, horizon T , counterfactual transition probability P τ , reward function R, initial state s 0 .
where note that the dynamics of the original states s are equivalent both under P + τ,t and P τ,t , however, under P + τ,t , we also keep track of the number of actions differing from the observed actions. Now, let π : S + × {0, . . . , T − 1} → A be a policy that deterministically decides about the counterfactual action a t that should have been taken if the process's enhanced state had been s + t = (s t , l t ), i.e., the counterfactual state at time t was s t after performing l t action changes. Then, given such a counterfactual policy π, we can compute the corresponding average counterfactual outcome as follows: where τ = {((s t , l t ), a t )} T −1 t=0 is a realization of the non-stationary counterfactual MDP M τ with a t = π((s t , l t ), t) and the expectation is taken over all the realizations induced by the transition probability P + τ and the policy π. Here, note that, if π((s t , 0), t) = a t for all t ∈ {0, . . . , T − 1}, thenō π (τ ) = o(τ ) matches the outcome of the observed realization.
Then, our goal is to find the optimal counterfactual policy π * τ that maximizes the counterfactual outcome subject to a constraint on the number of counterfactual actions that can differ from the observed ones, i.e., maximize πō π (τ ) subject to where a 1 , . . . , a T −1 is one realization of counterfactual actions and a 1 , . . . , a T −1 are the observed actions. The constraint guarantees that any counterfactual action sequence induced by the counterfactual transition probability P + τ and the counterfactual policy π can differ in at most k actions from the observed sequence. Finally, once we have found the optimal policy π * τ , we can sample a counterfactual realization of the process and the optimal counterfactual explanation using Algorithm 1.

Finding Optimal Counterfactual Explanations via Dynamic Programming
To solve the problem defined by Eq. 7, we break the problem into several smaller sub-problems. Here, the key idea is to compute the counterfactual policy values that lead to the optimal counterfactual outcome recursively by expanding the expectation in Eq. 6 for one time step. We start by computing the highest average cumulative reward h(s, r, c) that one could have achieved in the last r steps of the decision making process, starting from state S T −r = s, if at most c actions had been different to the observed ones in those last steps. For c > 0, we proceed recursively: and, for c = 0, we trivially have that: with s ∈ S, r ∈ {1, . . . , T }, c ∈ {1, . . . , k}, and h(s, 0, c) = 0 for all s and c. In Eq. 8, the first parameter of the outer maximization corresponds to the case where, at time T − r, the observed action a T −r is taken and the second parameter corresponds to the case where, instead of the observed action, the best alternative action is taken. By definition, we can easily conclude that h(s 0 , T, k) is the average counterfactual outcome of the optimal counterfactual policy π * τ , i.e., the objective value of the solution to the optimization problem defined by Eq. 7, and we can recover the optimal counterfactual policy π * τ by keeping track of the action chosen at each recursive step in Eq. 8 and 9. The overall procedure, summarized by Algorithm 2, uses dynamic programming-it first computes the values h(s, 1, c) for all s and c and then proceeds with the rest of computations in a bottom-up fashion-and has complexity O(n 2 mT k). Finally, we have the following proposition (proven by induction in Appendix B): The counterfactual policy π * τ returned by Algorithm 2 is the solution to the optimization problem defined by Eq. 7.

Experiments on Synthetic Data
In this section, we evaluate Algorithm 2 on realizations of a synthetic decision making process 2 . To this end, we first look into the average outcome improvement that could have been achieved if at most k actions had been different to the observed ones in every realization, as dictated by the optimal counterfactual policy. Then, we investigate to what extent the level of uncertainty of the decision making process influences the average counterfactual outcome achieved by the optimal counterfactual policy as well as the number of distinct counterfactual explanations it provides.
Experimental setup. We characterize the synthetic decision making process using an MDP with states S = {0, . . . , n − 1} and actions A = {0, . . . , m − 1}, where n = 20 and m = 10, and time horizon T = 20. For each state s and action a, we set the immediate reward equal to R(s, a) = s, i.e., the higher the state, the higher the reward. To set the values of the transition probabilities P (S t+1 | S t = s t , A t = a t ), we proceed as follows. First we pick one s ∈ S uniformly at random and we set a weight w s = 1 and then, for the remaining states s ∈ S \ s , we sample weights It is easy to check that, for each state-action pair (s t , a t ) at time t, s t+1 = s is most likely to be observed in the next timestep t + 1. The parameter α controls the level of uncertainty.
Then, we compute the optimal policy that maximizes the average outcome of the decision making process by solving Bellman's equation using dynamic programming [21] and use this policy to sample the (observed) realizations as follows. For each realization, we start from a random initial state s 0 ∈ S and, at each time step t, we pick the action indicated by the optimal policy with probability 0.95 and a different action uniformly at random with probability 0.05. This leads to action sequences that are slightly suboptimal in terms of average   : Influence that the level of uncertainty of a synthetic decision making process has on the average counterfactual outcomeō π * τ (τ ) achieved by the optimal counterfactual policy π * τ as well as on the number of distinct counterfactual explanations π * τ provides. In all panels, we set n = 20, m = 10 and d = 1,000 and, in each experiment, use 500 realizations from 10 different instances of the decision making process (50 realizations per instance), each with different w s . In panel (c), for each realization, we sample 100 counterfactual realizations and compute the average number of unique explanations across realizations. Shaded regions correspond to 95% confidence intervals. outcome 3 . Finally, to compute the counterfactual transition probability P τ,t for each observed realization τ , we follow the procedure described in Section 2 with d = 1,000 samples for each noise posterior distribution 4 .
Results. We first measure to what extent the counterfactual explanations provided by the optimal counterfactual policy π * τ would have improved the outcome of the decision making process. To this end, for each observed realization τ , we compute the relative difference between the average optimal counterfactual outcome and the observed outcome, i.e., (ō π * τ (τ ) − o(τ ))/o(τ ). Figure 2 summarizes the results for different values of k. We find that the relative difference between the average counterfactual outcome and the observed outcome is always positive, i.e., the sequence of actions specified by the counterfactual explanations would have led the process realization to a better outcome in expectation. However, this may not come as a surprise given that the counterfactual policy π * τ is optimal, as shown in Proposition 1. In this context, note that, in the worst case, the counterfactual policy π * τ would trivially repeat the observed action sequence, leading to an average counterfactual outcome equal to the observed outcome with probability 1. Moreover, we observe that, as the sequences of actions specified by the counterfactual explanations differ more from the observed actions (i.e., k increases), the improvement in terms of expected outcome increases.
Next, we investigate how the level of uncertainty α of the decision making process influences the average counterfactual outcome achieved by the optimal counterfactual policy π * τ as well as the number of distinct counterfactual explanations π * τ provides. Figure 3 summarizes the results, which reveal several interesting insights. As the level of uncertainty α increases, the average counterfactual outcome decreases, as shown in panel (a), however, the relative difference with respect to the observed outcome increases, as shown in panel (b). This suggest that, under high level of uncertainty, the counterfactual explanations may be more valuable to a decision maker who aims to improve her actions over time. However, in this context, we also find that, under high levels of uncertainty, the number of distinct counterfactual explanations increases rapidly with k. As a result, it may be preferable to use relatively low values of k to be able to effectively show the counterfactual explanations to a decision maker in practice.

Experiments on Real Data
In this section, we evaluate Algorithm 2 using real patient data from a series of cognitive behavioral therapy sessions. Similarly as in Section 5, we first measure the average outcome improvement that could have been achieved if at most k actions had been different to the observed ones in every therapy session, as dictated by the optimal counterfactual policy. Then, we look into individual therapy sessions and showcase how Algorithm 2, together with Algorithm 1, can be used to highlight specific patients and actions of interest for closer inspection 5 . Appendix D contains additional experiments benchmarking the optimal counterfactual policy against several baselines 6 .
Experimental setup. We use anonymized data from a clinical trial comparing the efficacy of hypnotherapy and cognitive behavioral therapy [25] for the treatment of patients with mild to moderate symptoms of major depression 7 . In our experiments, we use data from the 77 patients who received manualized cognitive behavioral therapy, which is one of the gold standards in depression treatment. Among these patients, we discard four of them because they attended less than 10 sessions. Each patient attended up to 20 weekly therapy sessions and, for each session, the dataset contains the theme of discussion, chosen by the therapist from a pre-defined set of themes (e.g., psychoeducation, behavioural activation, cognitive restructuring techniques). Additionally, a severity score is included, based on a standardized questionnaire [27], filled by the patient at each session, which assesses the severity of depressive symptoms. For more details about the severity score and the pre-defined set of discussion themes refer to Appendix C.
To derive the counterfactual transition probability for each patient, we start by creating an MDP with n = 5 states and m = 11 actions. Each state s ∈ S = {0, . . . , 4} corresponds to a severity score, where small numbers represent lower severity, and each action a ∈ A corresponds to a theme from the pre-defined list of themes that the therapists discussed during the sessions. Moreover, each realization of the MDP corresponds to the therapy sessions of a single patient ordered in chronological order and time horizon T ∈ {10, . . . , 20} is the number of therapy sessions per patient. Here, we denote the set of realizations for all patients as T .
In addition, to estimate the values of the transition probabilities, we proceed as follows. For every state-action pair (s i , a), we assume a n-dimensional Dirichlet(α i1 , . . . , α i1 ) prior on the probabilities p j | i,a = P (S t+1 = s j | S t = s i , A t = a), where α ij = 1 if j ∈ {i − 1, i, i + 1} and α ij = 0.01 otherwise. Then, if we observe c j transitions from state s i to each state s j after action a in the patients' therapy sessions T , we have that the posterior of the probabilities p j | i,a is a Dirichlet(α i1 + c 1 , . . . , α in + c n ). Finally, to estimate the value of the transition probabilities P (S t+1 | S t = s i , A t = a), we take the average of 100,000 samples from the posterior probability p j | i,a . This procedure sets the value of the transition probabilities proportionally to the number of times they appeared in the data, however, it ensures that all transition probability values are  Figure 4: Performance achieved by the optimal counterfactual policy π * τ in a series of real manualized cognitive behavioral therapy sessions T , where each realization τ ∈ T includes all the sessions of an individual patient sorted in chronological order. Panel (a) shows the distribution of the relative difference between the average counterfactual outcomeō π * τ (τ ) achieved by π * τ and the observed outcome o(τ ), i.e., (ō π * τ (τ ) − o(τ ))/o(τ ), for k = 3 (refer to Appendix E for additional results under other k values). Panels (b) and (c) show the average counterfactual outcomeō π * τ (τ ) achieved by π * τ and the average number of unique counterfactual explanations provided by each π * τ , averaged across patients, against the number of actions k differing from the observed ones. In panel (c), for each realization, the average number of unique counterfactual explanations provided by π * τ is estimated using 1,000 counterfactual realizations. In all panels, we set d = 1,000 and use data from 73 patients. Shaded regions correspond to 95% confidence intervals. non zero and transitions between adjacent severity levels are much more likely to happen. Moreover, we set the immediate reward for a pair of state and action (s, a) equal to R(s, a) = 5 − s ∈ {1, . . . , 5}, i.e., the lower the patient's severity level, the higher the reward. Here, if some state-action pair (s, a) is never observed in the data, we set its immediate reward to R(s, a) = −∞. This ensures that those state-action pairs never appear in a realization induced by the optimal counterfactual policy. Finally, to compute the counterfactual transition probability P τ,t for each realization τ ∈ T , we follow the procedure described in Section 2 with d = 1,000 samples for each noise posterior distribution.
Results. We first measure to what extent the counterfactual explanations provided by the optimal counterfactual policy π * τ would have improved each patient's severity of depressive symptoms over time. To this end, for each observed realization τ ∈ T corresponding to each patient, we compute the same quality metrics as in experiments on synthetic data in Section 5. Figure 4 summarizes the results. Panel (a) reveals that, for most patients, the improvement in terms of relative difference between the average optimal counterfactual outcomē o π * τ (τ ) and the observed outcome o(τ ) is rather modest. Moreover, panel (b) also shows that the absolute average optimal counterfactual outcomeō π * τ (τ ), averaged across patients, does not increase significantly even if one allows for more changes k in the sequence of observed actions. These findings suggest that, in retrospect, the choice of themes by most therapists in the sessions was almost optimal. That being said, for 20% of the patients, the average counterfactual outcome improves a ≥3 % over the observed outcome and, as we will later discuss, there exist individual counterfactual realizations in which the counterfactual outcome improves much more than 3%. In that context, it is also important to note that, as shown in panel (c), the growth in the number of unique counterfactual explanations with respect to k is weaker than the growth found in the experiments with synthetic data and, for k ≤ 4, the number of unique counterfactual explanations is smaller than 10. This latter finding suggests that, in practice, it may be possible to effectively show, or summarize, the optimal counterfactual explanations, a possibility that we investigate next.
We focus on a patient for whom the average counterfactual outcomeō π * τ (τ ) achieved by the optimal policy π * τ with k = 3 improves 9.5% over the observed outcome o(τ ). Then, using the policy π * τ , also with k = 3, and the counterfactual transition probability P τ , we sample multiple counterfactual explanations τ using Algorithm 1 and look at each counterfactual outcome o(τ ). Figure 5(a) summarizes the results, which show that, in most of these counterfactual realizations, the counterfactual outcome is greater than the observed outcome-if at most k actions had been different to the observed ones, as dictated by the optimal policy,  Figure 5: Insights on the counterfactual explanations provided by the optimal counterfactual policy π * τ for one real patient who received manualized cognitive behavioral therapy. Panel (a) shows the distribution of the counterfactual outcomes o(τ ) for the counterfactual realizations τ induced by π * τ and P τ . Panel (b) shows, for each time step, how frequently a counterfactual explanation changes the observed action as well as the observed severity level and the severity level in the counterfactual realization with the highest counterfactual outcome. Here, darker colors correspond to higher frequencies and higher severities. Panel (c) shows the action changes in the unique counterfactual explanations (green) provided by π * τ along with the mean of counterfactual outcomes (r) that each one achieves and how frequently (f) they appear across the counterfactual realizations. Here, the bottom row shows the observed actions that were changed by at least one of the counterfactual explanations. Refer to Appendix C for a definition of the actions (i.e., themes). In all panels, we set d = 1,000 and the results are estimated using 1,000 counterfactual realizations.
there is a high probability that the outcome would have improved. Next, we investigate to what extent there are specific time steps within the counterfactual realizations τ where π * τ is more likely to suggest an action change. Figure 5(b) shows that, for the patient under study, there are indeed time steps that are overrepresented in the optimal counterfactual explanations, namely t ∈ {10, 13, 16}. Moreover, the first of these time steps (t = 10) is when the patient had started worsening their depression after an earlier period in which they showed signs of recovery. Remarkably, we find that, in the counterfactual realization τ with the best counterfactual outcome, the worsening is mostly avoided. Finally, we look closer into the actual action changes suggested by the optimal counterfactual policy π * τ . Figure 5(c) summarizes the results, which reveal that π * τ recommends replacing some of the sessions on cognitive restructuring techniques (CRT) by behavioral activation (BHA) consistently across counterfactual realizations τ , particularly at the start of the worsening period. We discussed this recommendation with one of the researchers on clinical psychology who co-authored Fuhr et al. [25] and she told us that, from a clinical perspective, such recommendation is sensible since, whenever the severity of depressive symptoms is high, it is very challenging to apply CRT and instead it is quite common to use BHA. Appendix F contains additional insights about other patients in the dataset.

Conclusions
We have initiated the study of counterfactual explanations in decision making processes in which multiple, dependent actions are taken sequentially over time. Building on a characterization of sequential decision making processes using MDPs and the Gumbel-Max SCM, we have developed a polynomial time algorithm to find optimal counterfactual explanations. Using synthetic and real data from cognitive behavioral therapy, we have shown that the counterfactual explanations our algorithm finds can provide valuable insights to enhance sequential decision making under uncertainty.
Our work opens up many interesting avenues for future work, which may solve some of its limitations. For example, we have considered sequential decision making processes with discrete states and actions satisfying the Markov property. Although this setting may fit a plethora of real-world applications, it would be interesting to extend our work to continuous states and actions and/or semi-Markov processes. Moreover, we experimentally validated our method using a single real dataset. It would be valuable to evaluate counterfactual explanations generated by our algorithm using additional datasets from other (medical) applications. In this context, it would be worth to consider applications in which the true transition probabilities are due to a machine learning algorithm. Finally, the usefulness of the counterfactual explanations given by our algorithm crucially depends on the particular structural causal model we have focused on. To make such explanations more practical, it would be important to consider alternative structural causal models and carry out a user study in which the counterfactual explanations are shared with the human experts (e.g., therapists) who took the observed actions.

A Further related work
Our work also builds upon further related work on interpretable machine learning and counterfactual inference. In terms of interpretable machine learning, in addition to the work on counterfactual explanations for one-step decision making processes discussed in Section 1 [4][5][6], there is also a popular line of work focused on feature-based explanations [28][29][30]. Feature-based explanations highlight the importance each feature has on a particular prediction by a model, typically through local approximation. In this context, one particular type of feature-based explanations that have been relatively popular in explaining the action choices of reinforcement learning agents in simple gaming environments (e.g., Atari games) are those in the form of saliency maps [31][32][33]. In terms of counterfactual inference, the literature has a long history [34], however, it has primarily focused on estimating quantities related to the counterfactual distribution of interest such as, e.g., the conditional average treatment effect (CATE). More broadly, one may also draw connections between our work and automated planning [35,36] and off-policy evaluation in reinforcement learning [37][38][39].

B Proof of Proposition 1
Using induction, we will prove that the policy value π τ ((s, l), t ) set by Algorithm 2 is optimal for every s ∈ S, l ∈ {0, . . . , k}, t ∈ {0, . . . , T − 1} in the sense that following this policy maximizes the average cumulative reward h(s, r, c) that one could have achieved in the last r = T − t steps of the decision making process, starting from state S T −r = s, if at most c = k − l actions had been different to the observed ones in those last steps. Formally: subject to Recall that, a 1 , . . . , a T −1 are the observed actions and the counterfactual realizations a 1 , . . . , a T −1 are induced by the counterfactual transition probability P + τ and the policy π. We start by proving the induction basis. Assume that a realization has reached a state s + T −1 = (s, l) at time T − 1, one time step before the end of the process. If c = 0 (i.e., l = k), following Equation 9, the algorithm will choose the observed action π τ ((s, l), t ) = a T −1 and return an average cumulative reward h(s, 1, 0) = R(s, a T −1 ) + s ∈S P τ,T −1 (s | s, a T −1 ) h(s , 0, 0) = R(s, a T −1 ), where h(s , 0, 0) = 0 for all s ∈ S. Since no more action changes can be performed at this stage, this is the only feasible solution and therefore it is optimal.
If c > 0, since h(s , 0, c) = h(s , 0, c − 1) = 0 for all s ∈ S it is easy to verify that Equation 8 reduces to h(s, 1, c) = max a∈A R(s, a) and π τ ((s, l), t ) = argmax a∈A R(s, a) is obviously the optimal choice for the last time step. Now, we will prove that, for a counterfactual realization being at state s + t = (s, l) at a time step t < T − 1 (i.e., r = T − t , c = k − l), the maximum average cumulative reward h(s, r, c) given by Algorithm 2 is optimal, under the inductive hypothesis that the values of h(s , r , c ) already computed for r < r, c < c and all s ∈ S are optimal. Assume that the algorithm returns an average cumulative reward h(s, r, c) by choosing action π τ ((s, l), t ) = a while the optimal solution gives an average cumulative reward OP T s,r,c > h(s, r, c) by choosing an action a * = a. Here, by τ t = {((s t , l t ) , a t )} T −1 t=t we will denote realizations starting from time t with a t = π τ ((s t , l t ) , t) where π τ is the policy given by Algorithm 2 and we will use τ * t if the policy is optimal. Also, we will denote a possible next state at time t + 1, after performing action a, as (s , l ) where l = l + 1 if a = a t , l = l otherwise and, c = k − l . Similarly, after performing action a * , we will denote a possible next state as (s , l * ) where l * = l + 1 if a * = a t , l * = l otherwise and, c * = k − l * . Then, we get: where, in (a), we expand the expectation for one time step, in (b), we replace the average cumulative reward starting from time step t + 1 with h(s , r − 1, c ) and OP T s ,r−1,c * for the policy of Algorithm 2 and the optimal one respectively and, in (c), we replace OP T s ,r−1,c * with h(s , r − 1, c * ) due to the inductive hypothesis.
It is easy to see that, it can either be a * = a t with c * = c or a * ∈ A \ a t with c * = c − 1. If c = 0, following Equation 9, the algorithm will choose the observed action (i.e., a = a t ). This is the only feasible solution, since a * = a t would give c * = −1 and l * = k − c * = k + 1, which is not a valid state. Therefore, we get a = a * = a t , which is a a contradiction. If c > 0, because of the max operator in Equation 8, for the action a chosen by Algorithm 2, it necessarily holds that: which is clearly a contradiction.

C Additional Details about the Cognitive Behavioral Therapy Dataset
Each patient's severity of depression is measured using the standardized questionnaire PHQ-9 [27], which consists of 9 questions regarding the frequency of depressive symptoms (e.g., "Feeling tired or having little energy?") manifested over a period of two weeks. The patient has to answer each question by placing themselves on a scale ranging from 0 ("Not at all") to 3 ("Nearly every day"). The sum of those answers, ranging from 0 to 27, reflects the overall depression severity and it is usually discretized into five categories, corresponding to no depression (0 − 4), mild depression (5 − 9), moderate depression (10 − 14), moderately severe depression (15 − 19), severe depression (20 − 27) 8 . In our experiments, the states S = {0, . . . , 4} correspond to these five categories. Each session of cognitive behavioral therapy contains information about the topic of discussion between the patient and the therapist, among 24 pre-defined topics [25], with some of the topics having similar content. For example, there were 4 topics about "cognitive restructuring techniques" which, we observed that, some therapists merged and covered in 2 sessions. Here, we grouped the above topics into the following eleven broader themes: • STR -First session: Introduction, discussing expectations, getting to know each other, discussing the current symptoms / problems, current life situation.
• BIO -Biography: A look at biography, family and social frame of reference, school and professional development, emotional development, partnerships, important turning points or crises.
• PSE -Psychoeducation: Discuss symptoms of depression, recognize and understand connections between feelings, thoughts and behavior (depression triangle) based on a situation analysis from the current / last episode, causes of depression, develop a disease model, explain the treatment approach in relation to the model.
• BHA -Behavioural activation: Focus on behaviour, discuss the vicious circle (depression spiral), discuss list of pleasant activities, attention to life balance, if necessary improve the daily structure, recognizing and eliminating obstacles and problems.
• REV -Review: Review of the last sessions, collection of strategies learned so far, find suitable strategies for typical situations, draft a personal strategy plan, plan further steps.
• CRT -Cognitive restructuring techniques: Discuss influence of thoughts on feelings and actions, identify thought patterns, discuss influence of automatic thoughts / basic assumptions, check the validity of automatic thoughts.
• INR -Interactional competence: Self-assessment of your own self-confidence, discuss current interpersonal issues and derive goals, carry out role plays, transfer into everyday life.
• THP -Re-evaluation of thought patterns: Review, evaluate and rename basic assumptions, schemes and general plans.
• RLP -Relapse prevention: Explain the risk of relapse, discuss early warning symptoms, recognize risk situations, develop suitable strategies.
• END -Closing session: Finding a good end to the therapy, looking back on the last 5-6 months, parting ritual.
• EXT -Extra material: Sleep disorders, problem-solving skills, brooding module "When thinking doesn't help", discuss the influence of rumination on mood and impairments in everyday life, progressive muscle relaxation.
In our experiments, the actions A correspond to these broader themes. However, since the themes STR and END appeared only in the first (t = 0) and last (t = T − 1) time steps of each realization, we kept them fixed and we did not allow these themes to be used as action changes during the time steps t = {1, . . . , T − 2}. Average counterfactual outcome Random Greedy Noisy greedy Algorithm 2 Figure 6: Performance achieved by the optimal counterfactual policy π * τ given by Algorithm 2 and the baseline policies in the same series of cognitive behavioral therapy sessions T introduced in Section 6. The plot shows the average counterfactual outcome 1 T τ ∈Tō πτ (τ ) achieved by π * τ and the baseline policies, averaged over the set of observed realizations T , against the number of actions k differing from the observed ones. For each observed realization, the average counterfactual outcome is estimated using 1,000 counterfactual realizations. Here, we set d = 1,000 and use data from |T | = 73 patients. Shaded regions correspond to 95% confidence intervals.

D Performance Comparison with Baseline Policies
Experimental setup. In this section, we compare the average counterfactual outcome achieved by the optimal counterfactual policy, given by Algorithm 2, with that achieved by several baseline policies. To this end, we use the same experimental setup as in Section 6, however, instead of setting R(s, a) = −∞ for every unobserved pair (s, a), we set R(s, a) = 5 − s ∈ {1, . . . , 5}, similarly as for the observed pairs. This is because, otherwise, we observed that there were always realizations under the baselines policies for which the counterfactual outcome was −∞. In our experiments, we consider with the following baselines policies: • Random: At each time step t, the policy chooses the next action a * uniformly at random if l t < k and it chooses a * = a t otherwise.
• Greedy: At each time step t, being at state (s t , l t ), the policy chooses the next action a * greedily, i.e., if l t < k, then and, if l t = k, a * = a t .
• Noisy greedy: At each time step t, being at state (s t , l t ), it chooses the next action a * as follows. If l t < k, a * is given by Eq. 12 with probability 0.5 and a * = a t otherwise. If l t = k, a * = a t .
Results. Figure 6 shows the average counterfactual outcomes achieved by the optimal policy, as given by Algorithm 2, and the above baselines for different k values. The results show that, as expected, the optimal policy outperforms all the baselines across the entire range of k values and, moreover, the competitive advantage is greater for smaller k values. In addition, we also find that the performance of the random baseline policy drops significantly as k increases, since, as discussed in Section 6, the observed trajectories are close to optimal in retrospect and, differing from them causes the random policy to worsen the counterfactual outcome.  Figure 7: Empirical distribution of the relative difference between the average counterfactual outcomeō π * τ (τ ) achieved by the optimal counterfactual policy π * τ and the observed outcome o(τ ), i.e., P[(ō π * τ (τ ) − o(τ ))/o(τ )] for several values of k. Here, we consider the same series of cognitive behavioral therapy sessions T introduced in Section 6. In all panels, we set d = 1,000 and use data from 73 patients.

E Average counterfactual improvement for additional values of k
In this section, we use the same experimental setup as in Section 6 to measure to what extent the counterfactual explanations provided by the optimal counterfactual policy π * τ would have improved each patient's severity of depressive symptoms, under other values of k. To this end, for each observed realization τ ∈ T corresponding to each patient, we compute the relative difference between the average optimal counterfactual outcome and the observed outcome, i.e., (ō π * τ (τ ) − o(τ ))/o(τ ). Figure 7 summarizes the results, which show that, similarly as in Section 6, for most patients, the improvement in terms of relative difference between the average optimal counterfactual outcomeō π * τ (τ ) and the observed outcome o(τ ) is modest. That being said, we observe that, as the sequences of actions specified by the counterfactual explanations differ more from the observed actions (i.e., k increases), the improvement in terms of expected outcome presents a slight increase.

F Insights About Individual Patients
In this section, we provide insights about additional patients in the dataset. For each of these additional patients, we follow the same procedure in Section 6, i.e., we use Algorithm 1 with the policy π * τ , with k = 3, to sample multiple counterfactual explanations τ and look at the corresponding counterfactual outcomes o(τ ). Figure 8 summarizes the results, where each row corresponds to a different patient. The results reveal several interesting insights. For most of the patients, all of the counterfactual realizations lead to counterfactual outcomes greater or equal than the observed outcome (left column), however, the difference between the average counterfactual outcome and the observed outcome is relatively small. Notable exceptions are a few patients for whom there is a small probability that the counterfactual outcome is worse than the observed one (top row) as well as patients for whom the difference between the average counterfactual outcome and the observed outcome is high (bottom row). Additionally, we also find that the actual action changes suggested by the optimal counterfactual policies π * τ are typically concentrated in a few time steps across counterfactual realizations (right column), usually at the beginning or the end of the realizations. Action changes, severity vs. t Figure 8: Insights on the counterfactual explanations provided by the optimal counterfactual policy π * τ for five real patients who received manualized cognitive behavioral therapy. Each row corresponds to a different patient with an observed realization τ . The panels in the left column show the distribution of the counterfactual outcomes o(τ ) for the counterfactual realizations τ induced by π * τ and P τ . The panels in the right column show, for each time step, how frequently a counterfactual explanation changes the observed action as well as the observed severity level and the severity level in the counterfactual realization with the highest counterfactual outcome. Here, darker colors correspond to higher frequencies and higher severities. In all panels, we set d = 1,000, k = 3, and the results are estimated using 1,000 counterfactual realizations.