# Context-Specific Causal Discovery for Categorical Data Using Staged Trees

Manuele Leonelli
School of Human Sciences and Technology
IE University
Madrid, Spain
<EMAIL_ADDRESS>&Gherardo Varando
Image Processing Laboratory
Universitat de València
València, Spain
<EMAIL_ADDRESS>
###### Abstract

Causal discovery algorithms aims at untangling complex causal relationships
using observational data only. Here, we introduce new causal discovery
algorithms based on staged tree models, which can represent complex and non-
symmetric causal effects. To demonstrate the efficacy of our algorithms, we
introduce a new distance, inspired by the widely used structural
interventional distance, to quantify the closeness between two staged trees in
terms of their corresponding causal inference statements. A simulation study
highlights the efficacy of staged trees in uncovering complex, asymmetric
causal relationship from data and a real-world data application illustrates
their use in a practical causal analysis.

##  1 Introduction

One of the major tasks in all areas of science is to uncover causal
relationships between variables of interest. Since experimental data is in
many cases unavailable, this task often comes down to discover such
relationships using observational data only. This is usually referred to as
_causal discovery_.

One of the most common approaches in causal discovery, as well as in causal
analysis, is to represent causal relationships via _directed acyclic graphs_
(DAGs). In such graphs if there is an edge pointing from a variable to
another, then the former is a direct cause of the other (Pearl, 2009). The
literature on causal discovery using DAGs is now extensive (see Glymour et
al., 2019, for a rewiev). The most common approaches are the PC-algorithm
(Spirtes et al., 2000), greedy equivalence search (Chickering, 2002) and
functional causal models (Hoyer et al., 2008; Shimizu et al., 2006).

Although more efficient and scalable causal discovery algorithms are still
developed (e.g. Monti et al., 2020; Viinikka et al., 2020), most of the recent
literature has focused on continuous random variables only. Attention to
causal discovery for observational discrete data has been lately very limited
(see e.g. Cai et al., 2018; Cowell and Smith, 2014; Huang et al., 2018; Peters
et al., 2010, for exceptions). The aim of this paper is to discuss flexible
and powerful causal discovery algorithms for categorical data embedding
complex and non-symmetric variable relationships and to highlight their
efficacy.

Whilst most causal discovery is carried out via DAG models, here we consider
staged tree models (Collazo et al., 2018; Smith and Anderson, 2008) which,
differently to DAGs, can represent a wide array of non-symmetric causal
effects between discrete variables. Bayesian MAP structural learning
algorithms for this model class have been introduced (Collazo and Smith, 2016;
Cowell and Smith, 2014; Freeman and Smith, 2011) as well as score-based ones
(Leonelli and Varando, 2021; Silander and Leong, 2013). The latter have been
implemented in the open-source stagedtreees R package (Carli et al., 2020) and
are used henceforth.

Whilst there has been an extensive effort in determining causal effects
(Genewein et al., 2020; Görgen et al., 2015; Thwaites et al., 2010; Thwaites,
2013) and equivalence classes (Duarte and Solus, 2020; Görgen and Smith, 2018;
Görgen et al., 2018), only one causal discovery algorithm has been proposed by
Cowell and Smith (2014). Here we provide a first free implementation of the
dynamic programming algorithm of Cowell and Smith (2014) based on the
stagedtrees R package as well as a first extensive simulation study to assess
its effectiveness. We demonstrate that staged trees are extremely powerful in
discovering complex dependence structures and in general outperform DAGs for
discrete data.

We do this by introducing and studying a new measure of dissimilarity between
causal models tailored to the topology of staged trees and inspired by the
widely-used structural intervention distance (SID) of Peters and Bühlmann
(2015) which we henceforth call _context-specific intervention discrepancy_
(CID). Differently to SID which only accounts for symmetric causal
relationships, our CID can more generally consider the difference between two
causal models by accounting for complex, non-symmetric causal effects.

The code with the implemented methods and the simulation experiments is
available at
https://github.com/gherardovarando/context_specific_causal_discovery.

##  2 Staged trees

Let $[p]=\\{1,\dots,p\\}$ and $\bm{X}=(X_{i})_{i\in[p]}$ be categorical random
variables with joint mass function $P$ and sample space
$\mathbb{X}=\times_{i\in[p]}\mathbb{X}_{i}$. For $A\subset[p]$, we let
$\bm{X}_{A}=(X_{i})_{i\in A}$ and $\bm{x}_{A}=(x_{i})_{i\in A}$ where
$\bm{x}_{A}\in\mathbb{X}_{A}=\times_{i\in A}\mathbb{X}_{i}$. We also let
$X_{-A}=(X_{i})_{i\in[p]\setminus A}$.

Let $(V,E)$ be a directed, finite, rooted tree with vertex set $V$, root node
$v_{0}$ and edge set $E$. For each $v\in V$, let $E(v)=\\{(v,w)\in E\\}$ be
the set of edges emanating from $v$ and $\mathcal{C}$ be a set of labels.

######  Definition 1.

An $\bf X$-compatible staged tree is a triple $(V,E,\theta)$, where $(V,E)$ is
a rooted directed tree and:

  1. 1.

$V={v_{0}}\cup\bigcup_{i\in[p]}\mathbb{X}_{[i]}$;

  2. 2.

For all $v,w\in V$, $(v,w)\in E$ if and only if
$w=\bm{x}_{[i]}\in\mathbb{X}_{[i]}$ and $v=\bm{x}_{[i-1]}$, or $v=v_{0}$ and
$w=x_{1}$ for some $x_{1}\in\mathbb{X}_{1}$;

  3. 3.

$\eta:E\rightarrow\mathcal{L}=\mathcal{C}\times\cup_{i\in[p]}\mathbb{X}_{i}$
is a labelling of the edges such that $\eta(v,\bm{x}_{[i]})=(\kappa(v),x_{i})$
for some function $\kappa:V\to\mathcal{C}$.

If $\eta(E(v))=\eta(E(w))$ then $v$ and $w$ are said to be in the same
_stage_.

Therefore, the equivalence classes induced by $\eta(E(v))$ form a partition of
the internal vertices of the tree in _stages_.

$v_{0}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$X_{1}=0$$X_{1}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$

$w_{0}$$w_{1}$$w_{2}$$w_{3}$$w_{4}$$w_{5}$$w_{6}$$X_{1}=0$$X_{1}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$

Figure 1: An example of an $(X_{1},X_{2},X_{3})$-compatible (left) and an
$(X_{1},X_{3},X_{2})$-compatible (right) staged trees.

Definition 1 first constructs a rooted tree where each root-to-leaf path, or
equivalently each leaf, is associated to an element of the sample space
$\mathbb{X}$. Then a labeling of the edges of such a tree is defined where
labels are pairs with one element from a set $\mathcal{C}$ and the other from
the sample space $\mathbb{X}_{i}$ of the corresponding variable $X_{i}$ in the
tree. By construction, $\bf X$-compatible staged trees are such that two
vertices can be in the same stage if and only if they correspond to the same
sample space. Although staged trees can be more generally defined without
imposing this condition (see e.g. Collazo et al., 2018), henceforth, and as
common in practice, we focus on $\bf{X}$-compatible staged trees only (see
Leonelli, 2019, for an example of a non $\bf{X}$-compatible tree).

Figure 1 (left) reports an $(X_{1},X_{2},X_{3})$-compatible stratified staged
tree over three binary variables. The coloring given by the function $\kappa$
is shown in the vertices and each edge $(\cdot,(x_{1},\ldots,x_{i}))$ is
labeled with $X_{i}=x_{i}$. The edge labeling $\eta$ can be read from the
graph combining the text label and the color of the emanating vertex. The
staging of the staged tree in Figure 1 is given by the partition
$\\{v_{0}\\}$, $\\{v_{1},v_{2}\\}$, $\\{v_{3},v_{4}\\}$, $\\{v_{5}\\}$ and
$\\{v_{6}\\}$.

The parameter space associated to an $\bf X$-compatible staged tree
$T=(V,E,\eta)$ with labeling $\eta:E\rightarrow\mathcal{L}$ is defined as

$\Theta_{T}=\Big{\\{}\bm{\theta}\in\mathbb{R}^{\eta(E)}\;|\;\forall e\in
E,\theta_{\eta(e)}\in(0,1)\textnormal{ and }\forall v\sum_{e\in
E(v)}\theta_{\eta(e)}=1\Big{\\}}$ (1)

Let $\bm{l}_{T}$ denote the leaves of a staged tree $T$. Given a vertex $v\in
V$, there is a unique path in $T$ from the root $v_{0}$ to $v$, denoted as
$\lambda(v)$. For any path $\lambda$ in $T$, let $E(\lambda)=\\{e\in
E:e\in\lambda\\}$ denote the set of edges in the path $\lambda$.

######  Definition 2.

The _staged tree model_ $\mathcal{M}_{T}$ associated to the $\bf X$-compatible
staged tree $(V,E,\eta)$ is the image of the map

$\begin{array}[]{llll}\phi_{T}&:&\Theta_{T}&\to\Delta_{|\bm{l}_{T}|-1}\\\
&&\bm{\theta}&\mapsto\Big{(}\prod_{e\in
E(\lambda(l))}\theta_{\eta(e)}\Big{)}_{l\in\bm{l}_{T}}\end{array}$ (2)

An element of $\mathcal{M}(T)$ in Definition (2) identifies a joint
probability $P_{\bm{\theta}}$ with conditional distributions

$P_{\bm{\theta}}(X_{i}=x_{i}|X_{[i-1]}=\bm{x}_{[i-1]})=\theta_{\eta(\bm{x}_{[i-1]},\bm{x}_{[i]})}\quad\forall\bm{x}\in\mathbb{X}\text{
and }i\in[p].$

In the staged tree in Figure 1 the staging $\\{v_{1},v_{2}\\}$ implies that
the conditional distribution of $X_{2}$ given $X_{1}=0$, represented by the
edges emanating from $v_{1}$, is equal to the conditional distribution of
$X_{2}$ given $X_{1}=1$. It is apparent now that the vertex staging represents
conditional independence: the fact that $v_{1}$ and $v_{2}$ are in the same
stage implies that $X_{1}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X_{2}$. The staging given by the blue vertices
implies the context-specific independence (Boutilier et al., 1996)
$X_{3}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X_{2}|X_{1}=0$: the independence between
$X_{3}$ and $X_{2}$ holds only for one of the two levels of $X_{1}$. For such
a staged tree there is no equivalent DAG representation, since it embeds non-
symmetric conditional independences (Varando et al., 2021).

Two staged trees $T$ and $S$ are said to be statistically equivalent if they
induce the same models, that is $\mathcal{M}_{T}=\mathcal{M}_{S}$.

######  Definition 3.

Let $T=(V,E,\eta)$ be an $\bf X$-compatible staged tree. The _depth_ of a
vertex $v\in V$ equals the number of edges in $\lambda(v)$. The staged tree
$T$ is called _stratified_ if every two vertices in the same stage have the
same depth.

The staged tree in Figure 1 is stratified. Henceforth, and as common in
practice, we focus on stratified $\bf X$-compatible staged trees.

###  2.1 Causal models based on staged trees

We can define a causal staged tree model as a collection of interventional
distributions in the intuitive way: the joint distribution of $\bm{x}$ in the
intervened model is obtained by the product of the parameters in the
corresponding root-to-leaf path where we replace parameters corresponding to
intervened variables.

######  Definition 4.

A staged tree causal model induced by an $\mathbf{X}$-compatible staged tree
$T=(V,E,\eta)$ is the class of interventional distributions defined, for each
parameter vectors $\bm{\theta}\in\Theta_{T}$, as follows:

$\displaystyle
P_{\bm{\theta}}\left(\bm{X}=\bm{x}|\operatorname{do}(\bm{X}_{I}=\bm{z}_{I})\right)$
$\displaystyle=\prod_{i\not\in
I}\theta_{\eta(\bm{x}_{[i]},\bm{x}_{[i-1]})}\prod_{i\in I}\delta(x_{k},z_{k})$
(3)
$\displaystyle=\left\\{\begin{matrix}\frac{P_{\bm{\theta}}(\bm{X}=\bm{x})}{\prod_{i\in
I}P_{\bm{\theta}}(X_{i}=x_{i}|\bm{X}_{[i-1]}=\bm{x}_{[i-1]})}&\text{if
}\bm{x}_{I}=\bm{z}_{I}\\\ 0&\text{otherwise}\end{matrix}\right.$

In particular under the empty intervention we recover the observational
distribution $P_{\bm{\theta}}$.

We say that two staged trees $T$ and $S$ are causally equivalent if they
induce the same class of interventional distributions as in Definition 4.
Obviously two causally equivalent staged trees are also statistically
equivalent but not viceversa.

We have that, as for DAGs, intervening on a variable, only affects downstream
variables,

$\displaystyle P\left(X_{i}|\operatorname{do}(\bm{X}_{I}=\bm{z}_{I})\right)$
$\displaystyle=P\left(X_{i}|\operatorname{do}(\bm{X}_{I}=\bm{z}_{I\cap[i-1]})\right)$
(4)
$\displaystyle=\sum_{\bm{x}_{-I\cap[i-1]}}P\left(X_{i}|\bm{X}_{-I\cap[i-1]}=\bm{x}_{-I\cap[i-1]},\bm{X}_{I\cap[i-1]}=\bm{z}_{I\cap[i-1]}\right).$

Thus, in particular,

$P\left(X_{i}=x_{i}|\operatorname{do}(\bm{X}_{[i-1]}=\bm{z}_{[i-1]})\right)=P\left(X_{i}=x_{i}|\bm{X}_{[i-1]}=\bm{z}_{[i-1]}\right)=\theta_{\eta(\bm{x}_{[i-1]},\bm{x}_{[i]})}.$

##  3 Causal discovery algorithms

Causal discovery with staged tree models has been discussed by Cowell and
Smith (2014) and in Collazo et al. (2018).

In particular, Cowell and Smith (2014) used the order search introduced by
Silander and Leong (2013) together with an exhaustive search for the best
stages structures (the labelling $\eta$). In general, the dynamic programming
approach for order search of Silander and Leong (2013) can be coupled with
every score that can be decomposed across the depth of the staged tree and
with algorithms that estimate the stages structure of a variable given the set
of preceding ones. This method effectively finds the best scoring order,
reducing the needed computations compared to the naive procedure that would
score all possible ordering of the variables.

Collazo et al. (2018) correctly reports that causal discovery with staged
trees can be effectively carried only if coupled with a method to compute the
statistical equivalence class of a model. While an algebraic characterization
of equivalence classes (Görgen and Smith, 2018) and an algorithm to compute
them have been defined (Görgen et al., 2018), at the moment a practical
implementation is not available.

We thus propose to use the dynamic programming algorithm of Silander and Leong
(2013) together with two heuristic algorithms for stages structures discovery
under a given order: (1) a backward hill-climbing algorithm to minimize BIC
score and (2) a k-means clustering of square root probabilities (similar to
Silander and Leong (2013)). In both cases we use the BIC score as criteria for
selecting the best ordering of the variables (see Görgen et al., 2020, for
details).

We implement the dynamic programming order search in the framework of the
stagedtrees R package (Carli et al., 2020). In particular our implementation
can be coupled with all the staging recovery algorithms available in the
stagedtrees package. We refer to the Appendix for a detailed description of
the methods used.

##  4 Context-specific interventional discrepancy

Similarly to the structural interventional distance (SID) for DAGs (Peters and
Bühlmann, 2015), which counts the number of wrongly estimated interventional
distributions, we can define a context interventional discrepancy, with
respect to a reference staged tree ($T$). Such a discrepancy measures the
extent of the errors done in computing context-specific interventions using a
different staged tree ($S$).

######  Definition 5.

Let $T=(V,E,\eta)$ be an $\bm{X}$-compatible staged event tree and
$S=(W,F,\nu)$ an $\bm{X}_{\pi}$-compatible staged event tree, where $\pi$ is a
permutation of $[p]$. We define the context interventional discrepancy
$\operatorname{CID}(T,S)$ as,

$\operatorname{CID}(T,S)=\sum_{i\in[p]}\operatorname{CID}_{i}(T,S),$

where $\operatorname{CID}_{i}\left(T,S\right)$ is the proportion of contexts
$\bm{x}_{[i-1]}\in\mathbb{X}_{[i-1]}$ for which the interventional
distribution $P(X_{i}|\operatorname{do}(\bm{X}_{[i-1]}=\bm{x}_{[i-1]}))$ is
wrongly inferred by $S$ with respect to $T$. Precisely,
$P(X_{i}|\operatorname{do}(\bm{X}_{[i-1]}=\bm{x}_{[i-1]}))$ is wrongly
inferred by $S$ with respect to $T$ if there exists $P\in\mathcal{M}_{T}$ such
that

$P\left(X_{i}=x_{i}|\bm{X}_{[i-1]}=\bm{x}_{[i-1]}\right)\neq
P(X_{i}=x_{i}|\bm{X}_{I}\in\\{\bm{y}_{I}\in\mathbb{X}_{I}:\nu(\bm{y}_{K},w)=\nu(\bm{x}_{K},u)\\}),$

where $J=\\{j:j>i\text{ and }\pi(j)<\pi(i)\\}$, $I=\\{j:j<i\text{ and
}\pi(j)<\pi(i)\\}$ and $K=I\cup J$.

We detail in the Appendix the algorithm to compute the $\operatorname{CID}$
between two staged trees.

Intuitively, CID measures how much a different staged trees $S$ can be used to
compute interventional distributions of the type
$P(X_{i}|\operatorname{do}(\bm{X}_{[i-1]}=\bm{x}_{[i-1]}))$. We choose to
consider only marginal distributions under interventions on the preceding
variables in the true model $T$, instead of pairwise interventional
distributions of the form $P(X_{i}|\operatorname{do}(X_{j}=x_{j}))$ as in the
definition of SID (Peters and Bühlmann, 2015). This is because our goal is to
quantify the effect context-specific interventions. However, other measures of
discrepancy between staged tree causal models could not be alternatively
defined.

As an example of the computation of CID$(T,S)$, consider the two staged trees
in Figure 1, where the left tree is $T$ and the right one is $S$. We need to
determine which interventional distributions for the left staged tree are
wrongly inferred by the right one. For example, consider the intervention
$\operatorname{do}(X_{1}=0,X_{2}=1)$ and the distribution
$P(X_{3}|\operatorname{do}(X_{1}=0,X_{2}=1))$ for $P\in\mathcal{M}_{T}$. We
have that,

$P(X_{3}|\operatorname{do}(X_{1}=0,X_{2}=1))=P(X_{3}|X_{1}=0,X_{2}=1),$

and, $J=\emptyset,I=\\{1\\}$, thus, because $v_{3}$ and $v_{4}$ belong to the
same stage in $T$:

$P(X_{3}|X_{1}=0)=P(X_{3}|X_{1}=0,X_{2}=0)=P(X_{3}|X_{1}=0,X_{2}=1),$

and $P(X_{3}|\operatorname{do}(X_{1}=0,X_{2}=1))$ is then correctly inferred
by $S$. On the other hand, we have that
$P(X_{3}|\operatorname{do}(X_{1}=1,X_{2}=1))$ is wrongly inferred by $S$
because $P(X_{3},|X_{1}=1,X_{2}=1)\neq P(X_{3}|X_{1}=1)$ in general. To see
this, notice that

$P(X_{3}|X_{1}=1)=\sum_{x_{2}=0,1}P(X_{3}|X_{1}=1,X_{2}=x_{2})P(X_{2}=x_{2}|X_{1}=1),$

and $P(X_{3},|X_{1}=1,X_{2}=1)\neq P(X_{3}|X_{1}=1)$ if, for example, we
choose $P(X_{2}=x_{2}|X_{1}=1)=0.5$ and $P(X_{3}|X_{1}=1,X_{2}=0)\neq
P(X_{3}|X_{1}=1,X_{2}=1)$.

As another example, consider the interventional distribution
$P(X_{2}|\operatorname{do}(X_{1}=1))$. In this case we have $I=\\{1\\}$ and
$J=\\{3\\}$, and since vertices $v_{1},v_{2}$ are in the same stage (and so
are $\\{w_{5},w_{4}\\}$ and $\\{w_{3},w_{6}\\}$), we have that

$P(X_{2}|X_{1}=1)=P(X_{2})=P(X_{2}|X_{1}\in\\{0,1\\}),$

that is, $S$ correctly infers the interventional distribution
$P(X_{2}|\operatorname{do}(X_{1}=1))$ for every $P\in\mathcal{M}_{T}$.

Similar to SID, the context-specific intervention discrepancy is not
symmetric. The following proposition collects some properties of the newly
defined measure.

######  Proposition 1.

The following properties hold for $\operatorname{CID}$.

  1. i.

$\operatorname{CID}(T,S)=0$ for every pair of causally equivalent staged trees
$S,T$.

  2. ii.

If $\operatorname{CID}(T,S)=0$ then for every $P\in\mathcal{M}(T)$ we have
that $P(X_{i}|do(X_{[i-1]}=x_{[i-1]}))$ is correctly computed in $S$.

  3. iii.

If $\mathcal{M}(T)\subseteq\mathcal{M}(S)$ and $\pi$ is the identity, then
$\operatorname{CID}(T,S)=0$.

  4. iv.

If $\mathcal{M}_{T}$ is the fully independence model then $CID(T,S)=0$ for
every $X_{\pi}$-compatible staged tree $S$.

The proof of the above proposition can be found in the Appendix and follows
directly from Definition 5.

As described Collazo et al. (2018) and detailed in Varando et al. (2021),
every DAG model over a categorical random vector can be written as a staged
tree with variables ordered as one of the topological orders of the DAG. We
can thus always transform DAGs into statistically equivalent staged trees and
compute the CID between the two equivalent staged trees as explained above. In
order to compare CID and SID we perform a simulation study where we sample
uniformly DAGs over $5$ binary variables and compute their CID and SID. The
results are reported in the two dimensional density and scatter plot in Figure
2. We can see that there is a high correlation between the two measures thus
highlighting that CID is a sensible measure that could be used not only for
non-symmetric models, but also for symmetric ones based on DAGs.

Figure 2: CID and SID between randomly generated DAGs over $5$ binary
variables.

##  5 Simulation experiments

We perform a simulation study to evaluate the feasibility of the proposed
approach and to demonstrate its superiority with respect to the classical DAG
algorithms under the assumption that the true model is a staged tree.

We simulate data from randomly generated staged trees model with different
degrees of complexity: number of stages per variables ($k\in\\{2,3,4\\}$) and
number of possible levels per categorical variable ($l\in\\{2,3,4\\}$). We
consider models with $2,3,4$ and $5$ variables, and sample sizes ranging from
$100$ to $10000$ observations. For each parameters’ combination we perform
$100$ repetition of the experiment each time randomly shuffling the order of
the variables to eliminate any possible bias of the search heuristics.

Figure 3: Context interventional discrepancy (CID) and Kendall tau distance
(KD) between the estimated and true model.

We run the staged trees approach described in Section 3 using the backward
hill-climbing search (best_bhc) and the k-means heuristic (best_kmeans) with
the number of clusters fixed to $2$. Our two routines are compared to two
classical DAG learning algorithms such as tabu search (tabu) and max-min hill-
climbing (mmhc) (Tsamardinos et al., 2006), both implemented in the bnlearn R
package (Scutari, 2010).

Results are displayed in Figure 3, where the average CID and Kendall tau
distance between the variable orderings are plotted as a function of the
sample size $N$ for the case of $5$ variables with maximum $4$ levels.
Additional results and computational time of the different algorithms are
reported in the Appendix.

We can observe that, as expected, methods based on staged trees are able to
better recover the causal structure of the true model. Since the true models
are randomly generated staged trees we can expect that algorithms which search
for the best DAG are not able to recover the true relationships between the
variables. More specifically, we observe that the method based on the k-means
algorithm works very well when the number of stages per variable matches the
number of clusters ($k=2$), while, with respect to CID, its performance
degrades for $k=3,4$. The backward hill-climbing methods, instead, requires a
bigger sample size, but it is able to perform well even when the true model is
more complex.

Additionally, it is interesting to notice that both staged tree approaches are
able to recover well the causal order of the variables in all considered
scenarios. This is especially interesting for the k-means algorithm which
performs better than the backward hill-climbing method with respect to the
Kendall distance even when it is misspecified ($k>2$).

For completeness, in the Appendix, we report results for a simulation study
where the true model is a randomly generated DAG: as expected, in this case
DAG methods perform better than the proposed staged trees algorithms.

##  6 Real-world data

We next illustrate the use of staged trees to uncover causal relationships
using data from the 2014 survey “Aspects on everyday life" collected by ISTAT
(the Italian National Institute of Statistics) (ISTAT, 2014). The survey
collects information from the Italian population on a variety of aspects of
their daily lives. For the purpose of this analysis we consider five of the
many questions asked in the survey: do you practice sports regularly? (S =
yes/no); do you have friends you can count on? (F = yes/no); do you believe in
people? (B = yes/no); do you feel your neighborhood is degraded (D = yes/no);
do you watch TV? (T = yes/no). Instances with a missing answer were dropped,
resulting in 35870 answers to the survey.

For each of the possible variable orderings (120 in total), a staged tree is
learned with a backward hill-climbing algorithm optimizing the BIC score. In
Figure 4 the BIC of the best model for each order is reported. There are six
models having the same BIC and belonging to the same equivalence class.
Therefore in the data there is no information to favour one of these orders
over the others.

Figure 4: BICs of the best staged tree for every possible order of the ISTAT
dataset.

For each BIC-minimizing staged tree, we construct a DAG representation of the
tree using the methodology of Varando et al. (2021): an edge between two
variables represents the presence of (non-)symmetric dependence and the lack
of an edge is interpreted as a symmetric conditional independence statement
(as in standard DAG models). It can be shown that for all six orders the
resulting DAG is complete. Given these six DAGs we then construct a partial
DAG (PDAG) as follows: if in all six DAGs two variables are connected by an
edge according to the same direction, the PDAG has that directed edge;
otherwise, two variables are joined by an undirected edge. The resulting PDAG
is reported in Figure 5 (left).

SFBDT

SFBDT

Figure 5: PDAGs for the ISTAT dataset, constructed from staged trees (left)
and from standard DAG algorithms (right).

From the PDAG in Figure 5 (left) we can deduce that there is a dependence
between whether an individual practice sports (S), has friends to count on
(F), and believes in others (B), but no causal statements about them can be
made. The PDAG includes three edges into the vertex D. Therefore the data
supports the hypothesis that S, F and B have a causal effect on whether an
individual feels their neighborhood is degraded. Similarly, all previous
variables have a causal effect on whether an individual watches TV. Although
the PDAG tells us that such causal effects exist, it does not report the
specific form of non-symmetric dependence they entertain.

$v_{0}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$v_{7}$$v_{8}$$v_{9}$$v_{10}$$v_{11}$$v_{12}$$v_{13}$$v_{14}$$v_{15}$$v_{16}$$v_{17}$$v_{18}$$v_{19}$$v_{20}$$v_{21}$$v_{22}$$v_{23}$$v_{24}$$v_{25}$$v_{26}$$v_{27}$$v_{28}$$v_{29}$$v_{30}$S
= yesS = noF = yesF = noF = yesF = no B = yesB = noB = yesB = noB = yesB = noB
= yesB = noD = yesD = noD = yesD = noD = yesD = noD = yesD = noD = yesD = noD
= yesD = noD = yesD = noD = yesD = noT = yesT = noT = yesT = noT = yesT = noT
= yesT = noT = yesT = noT = yesT = noT = yesT = noT = yesT = noT = yesT = noT
= yesT = noT = yesT = noT = yesT = noT = yesT = noT = yesT = noT = yesT = noT
= yesT = no

Figure 6: Staged tree maximizing the BIC for the order (S,F,B,D,T).

In order to understand how the variables causally depend on each other, we can
refer to one of the BIC-minimizing staged trees. Figure 6 reports the learned
staged tree for the order (S,F,B,D,T). The staging over the variables D and T
(vertices $v_{7}$ to $v_{30}$) shows an highly asymmetric dependence structure
which could not be represented by a DAG model. For instance, the staging
$\\{v_{7},v_{11},v_{13}\\}$ and $\\{v_{8},v_{12},v_{14}\\}$ implies that, in
the context S = no, F has no causal effect on D. Furthermore, it implies that,
in the context F = yes, S has no causal effect on D. The fact that the
vertices $v_{15},v_{16},v_{19}$ and $v_{20}$ are in the same stage implies
that, in the context S = yes and B = yes, F and D have no causal effect on T.
In words, for individuals that practice sports and believe in others, whether
or not they have friends and live in degraded neighborhoods does not affect
whether or not they watch television. Due to space constraints we do not
discuss the full causal interpretation of the staged tree in Figure 6, but it
is apparent that the flexibility of the staging enables for the intuitive
representation of complex non-symmetric causal relationships learned from
data.

Figure 5 (right) further reports the PDAG learned using standard DAG learning
algorithms implemented in bnlearn (Scutari, 2010). Since DAGs cannot represent
non-symmetric types of dependence, the model embeds conditional independences
which are not supported by the data: the BIC of the best staged tree is
169754, whilst the one of the learned DAG is 169789. Furthermore, notice that
all the edges of the PDAG on the right are undirected and therefore standard
DAG algorithms do not actually find any causal relationship in this dataset.

##  7 Conclusions

We introduced and implemented causal discovery algorithms based on staged
trees which extend classic DAG models to account for complex, non-symmetric
causal relationships. In order to assess the effectiveness of staged trees in
causal reasoning, we defined a new discrepancy that measures the agreement
between the interventional distributions of two staged trees. Our simulation
experiments demonstrate that if the data is simulated from a staged tree
model, and therefore embeds non-symmetric relationships between variables,
staged trees outperform DAG models. Our real-world application further
highlights the need for non-symmetric models since staged trees, despite of
their complexity, outperform DAGs in terms of penalized model fit.

In general we demonstrated that staged tree models can be a valuable tool for
causal discovery and various directions for future work are possible. As
showed in the Appendix due to space constraints, the proposed methods obtain
good results also in the bivariate case which usually is impossible to tackle
with standard DAGs since the two possible directions belong to the same Markov
equivalence class. We are currently focusing on the derivation of theoretical
results about the identifiability of the causal order when non-symmetric
relationships between two variables are present.

Additional heuristics to learn the stages structures are currently being
developed, and similarly different strategies for learning variable ordering.
While the methods described in the present work obtain good results, they are
lacking in scalability and new heuristics are needed to tackle a larger number
of variables efficiently.

## Acknowledgments and Disclosure of Funding

Gherardo Varando’s work was partly funded by the European Research Council
(ERC) Synergy Grant “Understanding and Modelling the Earth System with Machine
Learning (USMILE)” under Grant Agreement No 855187.

## References

  * Boutilier et al. [1996] C. Boutilier, N. Friedman, M. Goldszmidt, and D. Koller.  Context-specific independence in Bayesian networks.  In _Proceedings of the 12th Conference on Uncertainty in Artificial Intelligence_ , pages 115–123, 1996.
  * Cai et al. [2018] R. Cai, J. Qiao, K. Zhang, Z. Zhang, and Z. Hao.  Causal discovery from discrete data using hidden compact representation.  _Advances in Neural Information Processing Systems_ , 2018:2666, 2018.
  * Carli et al. [2020] F. Carli, M. Leonelli, E. Riccomagno, and G. Varando.  The R package stagedtrees for structural learning of stratified staged trees.  _arXiv:2004.06459_ , 2020.
  * Chickering [2002] D. M. Chickering.  Optimal structure identification with greedy search.  _Journal of machine learning research_ , 3(Nov):507–554, 2002.
  * Collazo and Smith [2016] R. A. Collazo and J. Q. Smith.  A new family of non-local priors for chain event graph model selection.  _Bayesian Analysis_ , 11(4):1165–1201, 2016.
  * Collazo et al. [2018] R. A. Collazo, C. Görgen, and J. Q. Smith.  _Chain Event Graphs_.  Chapmann & Hall, 2018.
  * Cowell and Smith [2014] R. G. Cowell and J. Q. Smith.  Causal discovery through MAP selection of stratified chain event graphs.  _Electronic Journal of Statistics_ , 8(1):965–997, 2014.
  * Duarte and Solus [2020] E. Duarte and L. Solus.  Algebraic geometry of discrete interventional models.  _arXiv:2012.03593_ , 2020.
  * Freeman and Smith [2011] G. Freeman and J. Q. Smith.  Bayesian MAP model selection of chain event graphs.  _Journal of Multivariate Analysis_ , 102(7):1152–1165, 2011.
  * Genewein et al. [2020] T. Genewein, T. McGrath, G. Déletang, V. Mikulik, M. Martic, S. Legg, and P. A. Ortega.  Algorithms for causal reasoning in probability trees.  _arXiv:2010.12237_ , 2020.
  * Glymour et al. [2019] C. Glymour, K. Zhang, and P. Spirtes.  Review of causal discovery methods based on graphical models.  _Frontiers in Genetics_ , 10:524, 2019.
  * Görgen and Smith [2018] C. Görgen and J. Q. Smith.  Equivalence classes of staged trees.  _Bernoulli_ , 24(4A):2676–2692, 2018.
  * Görgen et al. [2015] C. Görgen, M. Leonelli, and J. Q. Smith.  A differential approach for staged trees.  In _European Conference on Symbolic and Quantitative Approaches to Reasoning and Uncertainty_ , pages 346–355, 2015.
  * Görgen et al. [2018] C. Görgen, A. Bigatti, E. Riccomagno, and J. Q. Smith.  Discovery of statistical equivalence classes using computer algebra.  _International Journal of Approximate Reasoning_ , 95:167–184, 2018.
  * Görgen et al. [2020] C. Görgen, M. Leonelli, and O. Marigliano.  The curved exponential family of a staged tree.  _arXiv:2010.15515_ , 2020.
  * Hoyer et al. [2008] P. Hoyer, D. Janzing, J. M. Mooij, J. Peters, and B. Schölkopf.  Nonlinear causal discovery with additive noise models.  _Advances in Neural Information Processing Systems_ , 21:689–696, 2008.
  * Huang et al. [2018] B. Huang, K. Zhang, Y. Lin, B. Schölkopf, and C. Glymour.  Generalized score functions for causal discovery.  In _Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pages 1551–1560, 2018.
  * ISTAT [2014] ISTAT.  Multiscopo istat – aspetti della vita quotidiana.  UniData - Bicocca Data Archive, Milano. Codice indagine SN147. Versione del file di dati 2.0, 2014.
  * Leonelli [2019] M. Leonelli.  Sensitivity analysis beyond linearity.  _International Journal of Approximate Reasoning_ , 113:106–118, 2019.
  * Leonelli and Varando [2021] M. Leonelli and G. Varando.  Structural learning of simple staged trees.  2021\.
  * Monti et al. [2020] R. P. Monti, K. Zhang, and A. Hyvärinen.  Causal discovery with general non-linear relationships using non-linear ica.  In _Uncertainty in Artificial Intelligence_ , pages 186–195. PMLR, 2020.
  * Pearl [2009] J. Pearl.  _Causality_.  Cambridge University Press, 2009.
  * Peters and Bühlmann [2015] J. Peters and P. Bühlmann.  Structural intervention distance for evaluating causal graphs.  _Neural Computation_ , 27(3):771–799, 2015.
  * Peters et al. [2010] J. Peters, D. Janzing, and B. Schölkopf.  Identifying cause and effect on discrete data using additive noise models.  In _Proceedings of the 13th International Conference on Artificial Intelligence and Statistics_ , pages 597–604. JMLR Workshop and Conference Proceedings, 2010.
  * Scutari [2010] M Scutari.  Learning Bayesian networks with the bnlearn R package.  _Journal of Statistical Software_ , 35(3):1–22, 2010.
  * Shimizu et al. [2006] S. Shimizu, P. O. Hoyer, A. Hyvärinen, A. Kerminen, and M. Jordan.  A linear non-Gaussian acyclic model for causal discovery.  _Journal of Machine Learning Research_ , 7(10), 2006.
  * Silander and Leong [2013] T. Silander and T. Y. Leong.  A dynamic programming algorithm for learning chain event graphs.  In _Proceedings of the International Conference on Discovery Science_ , pages 201–216, 2013.
  * Smith and Anderson [2008] J. Q. Smith and P. E. Anderson.  Conditional independence and chain event graphs.  _Artificial Intelligence_ , 172(1):42 – 68, 2008\.
  * Spirtes et al. [2000] P. Spirtes, C. N. Glymour, R. Scheines, and D. Heckerman.  _Causation, prediction, and search_.  MIT press, 2000.
  * Thwaites [2013] P. Thwaites.  Causal identifiability via chain event graphs.  _Artificial Intelligence_ , 195:291–315, 2013.
  * Thwaites et al. [2010] P. Thwaites, J.Q. Smith, and E. Riccomagno.  Causal analysis with chain event graphs.  _Artificial Intelligence_ , 174(12-13):889–909, 2010.
  * Tsamardinos et al. [2006] I. Tsamardinos, L. E. Brown, and C. F. Aliferis.  The max-min hill-climbing Bayesian network structure learning algorithm.  _Machine learning_ , 65(1):31–78, 2006.
  * Varando et al. [2021] G. Varando, F. Carli, and M. Leonelli.  Staged trees and asymmetry-labeled DAGs.  2021\.
  * Viinikka et al. [2020] J. Viinikka, A. Hyttinen, J. Pensar, and M. Koivisto.  Towards scalable Bayesian learning of causal DAGs.  In _Advances in Neural Information Processing Systems_ , volume 33, 2020.

##  Appendix A Structural learning of stratified staged trees

We report here the two algorithms to search the space of $\bm{X}$-compatible
staged trees available in the stagedtrees R package [Carli et al., 2020].
Moreover we review the dynamic programming approach to learn the best variable
ordering.

###  A.1 Learning the stage structure with a fixed order

The space of possible $\bm{X}$-compatible staged event tree is considerably
larger than the space of possible DAGs. Even if we fix the order of the
variables, exploring all possible combinations of stages structure becomes
rapidly infeasible [Collazo et al., 2018]. We thus use two of the possible
heuristic searches available in the stagedtrees package [Carli et al., 2020].

The backward hill-climbing method consist in starting from the saturated model
and, for each variable, iteratively try to join stages. At each step of the
algorithm all possible combination of two stages are tried and the best move
is chosen. The BIC score is used as guiding criterion; and since the log-
likelihood decomposes across the depth of the tree, the stages search can be
performed independently for each stratum.

The use of the k-means clustering of probabilities to learn staged event tree
was firstly introduced by Silander and Leong [2013] as a fast alternative to
backward hill-climbing algorithms. The default version implemented in the
stagedtrees package perform k-means clustering over the square-root of the
probabilities of a given variable given all the possible contexts.

Both algorithms operate over the stage structures of each variable $X_{i}$
independently of the other variable stages. Furthermore, the estimated stages
structure of a given variable is independent of the order of the preceding
variables and depends only on which variables precede $X_{i}$.

###  A.2 Learning an optimal variable ordering

The methods described in the previous section output an
$\bm{X}_{\pi}$-compatible staged tree for each possible ordering $\pi$ of the
variables. For a small number of variables, it is possible to simply enumerate
all possible $n!$ staged trees for all possible orders, and select the best
one(s) according to a chosen criterion (e.g. BIC).

Silander and Leong [2013] proposed a dynamic programming algorithm which still
obtains the global optimum, but with a substantial reduction in computational
complexity. The method can be coupled with every algorithm which operates
independently on every variable and using as guiding score every function
which can be decomposed across the variables of the model.

The search algorithm is based on the following assumptions:

  1. 1.

The estimated stages structures of a variable $X_{i}$ does not depend on the
order of the preceding variables. For example, the stage structure obtained
for variable $X_{3}$ (by one of the two methods in Section A.1) is the same
for the orderings $X_{1},X_{2},X_{3}$ and $X_{2},X_{1},X_{3}$.

  2. 2.

The estimated stages structure of a variable is independent of the estimated
stages structures for the remaining variables.

  3. 3.

The selection criterion to be optimized decomposes across the variables.

We refer to Silander and Leong [2013] and Cowell and Smith [2014] for the
specifics of the dynamic programming approach.

##  Appendix B Proof of Proposition 1

If two staged tree $T,S$ are causally equivalent then for every
$P\in\mathcal{M}_{T}=\mathcal{M}_{S}$ we have

$P(X_{i}|\operatorname{do}(\bm{X}_{[i-1]}=\bm{x}_{[i-1]}))=P(X_{i}|\operatorname{do}(\bm{X}_{I}=\bm{x}_{I})),$

thus,

$\displaystyle P(X_{i}=x_{i}|\bm{X}_{[i-1]}=\bm{x}_{[i-1]})$
$\displaystyle=P(X_{i}=x_{i}|\operatorname{do}(\bm{X}_{I}=\bm{x}_{I}))$
$\displaystyle=P(X_{i}=x_{i}|\bm{X}_{I}\in\\{\bm{y}_{I}\in\mathbb{X}_{I}:\nu(\bm{y}_{K},w)=\nu(\bm{x}_{K},u)\\})$

which proves point (i). Point (ii) follows directly from the definition of
CID.

To prove point (iii), observe that, since $\pi$ is the identity both
$T=(V,E,\eta)$ and $S=(V,E,\nu)$ are $\bm{X}$-compatible staged event trees
and thus $\mathcal{M}_{T}\subseteq\mathcal{M}_{S}$ implies that
$\nu(u,v)=\nu(u^{\prime},v^{\prime})\Rightarrow\eta(u,v)=\eta(u^{\prime},v^{\prime})$
(the stage structure of $T$ is coarser than the one of $S$). Since
$P\in\mathcal{M}_{T}$, we have,

$P(X_{i}=x_{i}|\bm{X}_{[i-1]}=\bm{x}_{[i-1]})=P(X_{i}=x_{i}|\bm{X}_{[i-1]}\in\\{\bm{y}_{[i-1]}\in\mathbb{X}_{[i-1]}:\nu(\bm{y}_{[i-1]},w)=\nu(\bm{x}_{[i-1]},u)\\}).$

Finally, point (iv) follows from points (i) and (iii) by observing that if $T$
is the completely independent model then $T$ is causally equivalent to any
completely independent $\bm{X}_{\pi}$-compatible staged tree, for any
permutation $\pi$.

##  Appendix C Algorithm for CID

The algorithm to compute the context-specific interventional discrepancy is
reported in Algorithm 1.

In the pseudo-code for the algorithm we use the following notation, for a
staged event tree $S=(W,F,\nu)$ we denote with $\sim_{S}$ the equivalence
relation over $V$ defined by $u\sim_{S}v$ if and only if
$\nu(E(v))=\nu(E(v))$. Thus the equivalence classes with respect to the above-
defined relation are the stages of $S$.

Algorithm 1 Compute context interventional discrepancy

0: $T=(V,E,\eta)$ an $\bm{X}$-compatible staged event tree and $S=(W,F,\nu)$
an $\bm{X}_{\pi}$-compatible staged event tree.

0: The CID between $T$ and $S$.

initialize $\texttt{CID}=0$

for $i=1$ to $p$ do

$k=\pi^{-1}(i)$ # the position of variable $X_{i}$ in staged tree $S$

$I=\\{j:j<i\And\pi^{-1}(j)<k\\}$

$J=\\{i:j>i\And\pi^{-1}(j)<k\\}$

$\mathcal{A}=\mathbb{X}_{\pi([k-1])}/\sim_{S}$ # the stages of $S$ at depth
$i$

$\texttt{wrong}=\emptyset$

for $A\in\mathcal{A}$ do

$B_{A}=\\{\bm{x}_{[i-1]}\in\mathbb{X}_{[i-1]}:\bm{x}_{I}=\bm{y}_{I}\text{ for
some }\bm{y}\in A\\}$

if $|\eta(E(B_{A}))|>1$ then

$\texttt{wrong}=\texttt{wrong}\cup B_{A}$

end if

$\texttt{CID}=\texttt{CID}+\frac{|\texttt{wromg}|}{|\mathbb{X}_{[i-1]}|}$

end for

end for

We prove now that Algorithm 1 correctly computes the CID.

######  Proposition 2.

For each pair of $X$ and $X_{\pi}$ compatible staged event trees, $T$ and $S$,
the output of Algorithm 1 is equal to $CID(T,S)$.

###### Proof.

We need only to prove that the procedure in the first loop effectively
identifies wrongly inferred interventional distribution of the type
$P(X_{i}|\operatorname{do}(\bm{X}_{[i-1]}=\bm{x}_{[i-1]}))$. That is, at the
end of $i$-th iteration of the first for loop, the variable wrong contains the
set of contexts $\bm{x}_{[i-1]}\in\mathbb{X}_{[i-1]}$ such that the
corresponding interventional distribution for $X_{i}$ is wrongly inferred by
$S$ with respect to $T$. If $\bm{x}_{[i-1]}\in\texttt{wrong}$ then, by
construction, there exists $A\in\mathcal{A}$ such that $\bm{x}_{[i-1]}\in
B_{A}$, and there exists $\bm{x}^{\prime}_{[i-1]}\in B_{A}$ with
$\eta(E(\bm{x}_{[i-1]}))\neq\eta(E(\bm{x}^{\prime}_{[i-1]}))$. Thus, there
exists $P\in\mathcal{M}_{T}$ such that
$P(X_{i}|\bm{X}_{[t-1]}=\bm{x}_{[i-1]})\neq
P(X_{i}|\bm{X}_{[t-1]}=\bm{x}^{\prime}_{[i-1]})$ and moreover

$P(X_{i}|\bm{X}_{[t-1]}=\bm{x}_{[i-1]})\neq P(X_{i}|X_{i}\in B_{A}).$

Viceversa, if $\bm{x}_{[i-1]}\not\in\texttt{wrong}$, for every
$P\in\mathcal{M}_{T}$ we have that, $\eta(E(B_{A}))$ is a singleton and thus,

$P(X_{i}|\bm{X}_{[t-1]}=\bm{x}_{[i-1]})=P(X_{i}|X_{i}\in B_{A}),$

for every $A\in\mathcal{A}$ such that $\bm{x}_{[i-1]}\in B_{A}$. And thus $S$
correctly infers $P(X_{i}|\bm{X}_{[i-1]}=\bm{x}_{[i-1]})$. ∎

##  Appendix D Additional simulation results

For completeness we report here additional simulation results, which for space
constraints could not fit in the main paper. We report the CID and Kendall tau
distance between the true model and the estimated ones. The Kendall tau
distance is the computed between the true causal order and the estimated order
with the implementation in the PerMallows R package [permallows]. The full
simulations can be replicated following the instructions available at
https://github.com/gherardovarando/context_specific_causal_discovery.

Figure 7 reports the computational times of the algorithms. As expected,
algorithms for DAGs are faster since the searched model space is much smaller.
The k-means algorithm for staged trees is comparable to those for DAGs in
terms of speed and its complexity does not seem to exponentially increase as
in the case of the backward hill-climbing.

Figure 7: Computational time in seconds for different algorithms. Figure 8:
Simulation results for data sampled from randomly generated DAGs ($p=4$).

In Figure 8 we can see results of simulation where the true model is a
randomly generated DAG. As expected the tabu search over DAGs and the Max-Min
Hill-Climbing algorithm (mmhc) [Scutari, 2010] obtain better results then the
staged trees methods. Still, we noticed in the paper that the causal order is
retrieved poorly by DAGs algorithms when the true model is a staged event
tree. We can thus hypothesise that the asymmetric conditional independence
statements present in a random staged event tree model are helpful to retrieve
the causal structure.

Interestingly, when the true model is a randomly generated staged event tree,
even in the bivariate case, we observe that the proposed staged trees methods
outperform classical DAG methods 9. This suggests that the asymmetrical stages
structure in a staged event tree can be sufficient to have identifiability of
the causal order. In the next section we report a concrete example where such
identifiability is explicit.

Figure 9: Simualtion results for the bivariate case ($p=2$). Data sampled from
randomly generated staged event trees (with $k$ stages per variable).

##  Appendix E An example of identifiability of the causal order

As an instructive example we report here a bivariate staged tree model
compatible with $(X_{1},X_{2})$ where the causal order can be identified by
choosing the simpler model.

Consider the two staged tree depicted in Figure 10. Let $T$ be the staged tree
on the right and $S$ the one on the left.

It is easy to see that if the data are generated from a joint probability
distribution $P\in\mathcal{M}_{T}$, $S$ is the only $(X_{1},X_{2})$-compatible
staged tree such that $P\in\mathcal{M}_{S}$. Indeed, since $\mathcal{M}_{S}$
is the saturated model we have that $\mathcal{M}_{S}$ is equal to the entire
probability simplex.

But we can see that $\mathcal{M}_{T}\subsetneq\mathcal{M}_{S}$, since
$P(X_{1}|X_{2}=2)=P(X_{1}|X_{2}=3)$ in $\mathcal{M}_{T}$, and thus the causal
order is here identified by choosing the simplest model which describes the
data generating process. Even if this is a very simple example, it is
instructive to see that non-symmetrical conditional independence statements
can be leveraged by staged event trees to discover causal structure in
categorical data.

$v_{0}$$v_{1}$$v_{2}$$X_{1}=0$$X_{1}=1$$X_{2}=1$$\quad
X_{2}=2$$X_{2}=3$$X_{2}=1$$\quad
X_{2}=2$$X_{2}=3$$w_{0}$$w_{1}$$w_{2}$$w_{3}$$X_{2}=1$$X_{2}=2$$X_{2}=3$$X_{1}=0$$X_{1}=1$$X_{1}=0$$X_{1}=1$$X_{1}=0$$X_{1}=1$
Figure 10: An example of an $(X_{1},X_{2})$-compatible (left) and an
$(X_{2},X_{1})$-compatible (right) staged trees.


# Detecting Post-Stroke Aphasia Via Brain Responses to Speech in a Deep
Learning Framework

Pieter De Clercq11, Corentin Puffay11,2, Jill Kries1, Hugo Van Hamme2,
Maaike Vandermosten1, Tom Francart1 and Jonas Vanthornhout1 1Shared first
authors
1KU Leuven, dept. Neurosciences, ExpORL, Leuven, Belgium
2KU Leuven, dept. Electrical engineering (ESAT), PSI, Leuven, Belgium
Correspondence<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Aphasia, a language disorder primarily caused by a stroke, is traditionally
diagnosed using behavioral language tests. However, these tests are time-
consuming, require manual interpretation by trained clinicians, suffer from
low ecological validity, and diagnosis can be biased by comorbid motor and
cognitive problems present in aphasia. In this study, we introduce an
automated screening tool for speech processing impairments in aphasia that
relies on time-locked brain responses to speech, known as neural tracking,
within a deep learning framework. We modeled electroencephalography (EEG)
responses to acoustic, segmentation, and linguistic speech representations of
a story using convolutional neural networks trained on a large sample of
healthy participants, serving as a model for intact neural tracking of speech.
Subsequently, we evaluated our models on an independent sample comprising 26
individuals with aphasia (IWA) and 22 healthy controls. Our results reveal
decreased tracking of all speech representations in IWA. Utilizing a support
vector machine classifier with neural tracking measures as input, we
demonstrate high accuracy in aphasia detection at the individual level
(85.42%) in a time-efficient manner (requiring 9 minutes of EEG data). Given
its high robustness, time efficiency, and generalizability to unseen data, our
approach holds significant promise for clinical applications.

##  I Introduction

The prevalence of aphasia, a language disorder primarily resulting from a
stroke, is substantial, impacting over one-third of global stroke patients
[1]. Accurate and timely diagnosis of aphasia typically performed using
behavioral language tests, is crucial for initiating effective therapy.
However, language tests used in the clinic do not measure natural, daily-life
language skills, as they consist of artificial tasks that target isolated
language components (i.e., phonetic, semantic, or syntactic processes) [2].
Furthermore, these tests are time-consuming in administration and post-hoc
interpretation by a trained clinician and can be biased by comorbid motor and
cognitive problems present in over 80% of individuals with aphasia (IWA) [3,
4].

To address these limitations, electroencephalography (EEG) studies have
explored biomarkers of impaired language functioning in aphasia. EEG-based
neural tracking of natural speech measures time-locked brain responses to
acoustic (e.g., the speech envelope) and linguistic (providing the semantic
value carried by a word or a phoneme) [5] speech representations. Several
methods are used to relate the neural signals to speech, for example using
backward models that reconstruct the speech features from the neural signals
[6] or models that attempt to discriminate time-aligned from misaligned speech
segments (i.e., the match-mismatch paradigm) using the neural signals [7].
Recent findings using linear models showed decreased neural tracking of the
speech envelope and segmentation features (i.e., phoneme and word onsets) in
IWA compared to healthy controls, but no group differences were observed for
linguistic speech tracking [8]. Focusing on the temporal envelope of speech, a
follow-up study using nonlinear mutual information analyses demonstrated that
primarily low-frequency envelope tracking (delta and theta band) is impaired
in IWA [9], which carries cues for detecting and segmenting lexical units of
speech (i.e., syllables, words, phrases and sentences) and prosody [10]. Using
neural envelope tracking measures as input, classification algorithms reliably
detected aphasia at the individual level [9]. Together, these results show
promise for neural tracking as an automated assessment tool for natural speech
abilities in aphasia.

A critical drawback in previous aphasia studies is that the applied techniques
(i.e., linear models/mutual information) require participant- and stimulus-
specific training data. Here, we test the potential of more powerful deep
neural networks to generalize to unseen data and capture speech processing
impairments in IWA, aiming to build a more robust and generalizable screening
tool for aphasia. Specifically, we pre-train convolutional neural networks
(CNNs) on acoustic (envelope), segmentation (phoneme/word onsets), and
linguistic speech representations in a large independent sample of healthy
individuals with intact speech processing. Subsequently, we evaluate the pre-
trained CNNs on individual data from IWA and healthy age-matched controls. All
models’ performance for each speech feature is computed and provided as input
to a nonlinear support vector machine (SVM) classifier trained to detect
aphasia at the individual level. Furthermore, we investigate the contribution
of each feature and the minimally required recording length for accurate
aphasia detection.

##  II Methods

###  II-A Participants

The presented sample is identical to previous work in our lab [9], comprising
26 IWA (seven female participants, mean age=72, std=15 y/o) in the chronic
phase ($>$ 6 months after stroke) and 22 neurologically healthy controls
(seven female participants, mean age=72, std=7 y/o), all Flemish-native
speakers. IWA suffered from a left-hemispheric or bilateral stroke and were
diagnosed with aphasia in the acute stage post-stroke. IWA and healthy
controls both performed standardized clinical tests for aphasia at the time of
study participation, with IWA displaying worse performance [9]. Although 27%
of IWA did not score below threshold to be categorized as aphasia on either of
these tests, these patients had an extended language deficit documentation and
followed speech-language therapy at the time of study participation. Hence,
they were still labelled as aphasia. For more details on recruitment strategy
and the aphasia sample, we refer to [9]. The study was approved by the ethical
committee UZ/KU Leuven (S60007), and all participants gave written consent
before participation.

During the EEG experiment, participants were instructed to listen to a
25-minute-long story, ”De Wilde Zwanen”, written by Christian Andersen and
narrated by a female Flemish-native speaker. For each subject, we used 20
minutes of the entire recording, as we experienced data loss for 1 participant
in the last five minutes of recording. EEG data was recorded using a
64-channel Active-Two EEG system (BioSemi, Amsterdam, Netherlands) at 8192 Hz
sampling rate.

###  II-B Signal processing and feature extraction

####  II-B1 EEG data processing

We performed pre-processing steps on the EEG data, including downsampling to
512 Hz, independent component analysis to remove artifacts, common average-
referencing and filtering in frequency bands of interest: delta (0.5-4 Hz),
theta (4-8 Hz), alpha (8-12 Hz), beta (12-25 Hz) and a broad (0.5-32 Hz) band.
We used a Least Squares filter of order 2000 with 10% transition bands
(transition of frequencies 10% above the lowpass filter and 10% below the
highpass) and compensation for the group delay. Furthermore, the EEG signals
were z-score normalized and downsampled to 64 Hz.

####  II-B2 Speech features

We extracted the speech features from the stimulus to investigate neural
tracking at the acoustic (5 features), segmentation (2 features), and
linguistic (4 features) level. All features were downsampled to 64 Hz to be
synchronized with the EEG signal. For detailed information about the
implementation of features, we refer to our recent publications using the same
feature sets [11, 8, 9].

  * •

Acoustic features: We extracted the temporal envelope of speech in the delta,
theta, alpha, beta and broad band (reportedly involved in different stages of
speech processing [10]), following the procedure used in our prior aphasia
work [9]. The envelopes were filtered to the frequencies of interest using the
same Least Squares filter as applied on the EEG.

  * •

Segmentation features: The onsets of phonemes and words, were extracted from
the stimulus. The resulting features were one-dimensional vectors identifying
the presence of a word (or a phoneme) onset with an impulse of unit magnitude,
and zero otherwise.

  * •

Linguistic features:  The linguistic features are encoded identically to the
segmentation features, except for the unit magnitude that is replaced with a
value reflecting the information carried by the word or the phoneme. At the
phoneme level, we computed: (1) Phoneme surprisal: reflects the probability of
each phoneme given the preceding phonemes in the same word, and (2) Cohort
entropy: reflects the degree of competition among possible words created from
the phoneme of interest. At the word level, we computed (3) Word frequency:
reflects how frequently a word occurs in the language and, (4) Word surprisal:
reflects the probability of each word given the preceding words in the same
sentence.

###  II-C The match-mismatch task

We trained all our models on the supervised match-mismatch task [7]. As
depicted in Figure 1A, our models are trained to associate the EEG segment
with the matched speech segment among two presented labeled speech segments.
The matched speech segment is synchronized with the EEG, while the mismatched
speech segment occurs 1 s after the end of the matched segment, following
suggestions in prior work [12]. We fixed the segment length to 5 s for the
speech envelope and phoneme-based features. For word-based features, we chose
10 s segments to provide enough context to the models as suggested in [11].

Figure 1: A. Match-mismatch task. Finding which of the matched (yellow) or
mismatched (black) segments is synchronized with the EEG segment. B. Model
architectures. All models contain an EEG segment, a matched, and a mismatched
segment of speech. They are fed to a CNN block, compared with cosine
similarity, and passed through a last dense layer to select the matched
segment. The single-feature architecture contains only the top stream,
indicated with the red rectangle, while the dual-feature architecture includes
the entire scheme.

###  II-D Model architectures

To evaluate the neural tracking of speech features, we used two CNN-based
architectures depicted in Figure 1B: (1) the single-feature model [13, 14],
and (2) the dual-feature model [14, 11]. We used the single-feature
architecture for the broadband envelope and the different frequency bands
following up on [13]. For the segmentation and linguistic features, we used
the dual-feature model as proposed in [11]. The single-feature model consists
of two streams: one for the EEG segment, and one for the matched and
mismatched speech segments. The EEG stream contains one spatial convolution
followed by three temporal convolutions. The speech stream contains only the
three temporal convolutions. The embedded representation of the EEG is
compared using the cosine similarity with the embedded representation of the
matched and mismatched speech segments, and the labels (i.e., match or
mismatch) are predicted using a dense layer with a sigmoid activation
function. The dual-feature model is the duplicated single-feature model’s
architecture, using two speech features instead of one, to select the matched
segment. Both architectures are trained with the Adam optimizer and binary
cross-entropy as the loss function. For more details, we refer to [13, 14,
11]. We report in Table I the features used for the given architectures.

TABLE I: Architectures used for each speech feature. Speech feature | Architecture
---|---
Envelopes (N=5) | single-feature
2 x word onset | dual-feature
2 x phoneme onset | dual-feature
Word frequency + word onset | dual-feature
Word surprisal + word onset | dual-feature
Phoneme surprisal + phoneme onset | dual-feature
Cohort entropy + phoneme onset | dual-feature

###  II-E Model pre-training and evaluation on the match-mismatch task

Each model we used was pre-trained on the match-mismatch task with a dataset
consisting of 60 young native Flemish subjects listening to various speech
materials involving different speakers, as used previously [11]. We trained,
validated, and tested each model using an 80/10/10 split for each subject,
following the method of our previous publications [13, 15]. Subsequently, we
evaluated our models on the aphasia and age-matched control participants. For
each model, the evaluation metric is the average match-mismatch accuracy
across all segments per participant. This resulted in 11 accuracy values per
participant, one for each speech feature. We performed group comparisons on
accuracy values using non-parametric Wilcoxon rank-sum tests, and p-values
were adjusted for multiple comparisons employing the false discovery rate
(FDR) correction with the Benjamini-Hochberg procedure [16].

###  II-F Aphasia classification

The match-mismatch accuracies of all 11 models for each participant were used
to classify the aphasia status of participants. In addition, age (which is
matched between groups) was included as a feature, as it can affect neural
tracking outcomes [17]. We used a nonlinear radial basis function kernel SVM
and performed a nested cross-validation approach. In the inner cross-
validation, the C-hyperparameter and pruning were optimized (accuracy-based)
and tested in a validation set using 5-fold cross-validation. Predictions were
made on the test set in the outer loop using leave-one-subject-out cross-
validation. We computed the receiver operating characteristic (ROC) curve and
calculated the area under the curve (AUC), and further reported the accuracy,
F1-score, sensitivity and specificity of the classifier.

Model interpretation was performed using SHAP (SHapley Additive exPlanations)
analysis [18]. SHAP values quantify the impact of individual features on the
SVM’s classification, with positive values denoting features pushing
predictions towards the positive class (aphasia), vice-versa for the negative
class (healthy control), and values around 0 reflecting a limited feature
impact.

Finally, we determined the required EEG recording length for accurate aphasia
detection. We gradually increased the amount of EEG and speech data, ranging
from 1 min to 20 min, with a step size of 1 min. For each duration and
participant, we calculated the 11 match-mismatch accuracy values. We then
repeated the SVM procedure for aphasia classification and computed the
classification accuracy obtained for each recording length.

##  III Results

IWA exhibited statistically significant lower match-mismatch accuracies for
all features, as depicted in Figure 2: broadband ($W=78$, $p<.001$), delta
($W=160$, $p=0.009$), theta ($W=75$, $p<.001$), alpha ($W=64$, $p<.001$), beta
($W=150$, $p=0.005$), phoneme onsets ($W=90$, $p<.001$), word onsets ($W=117$,
$p<.001$), cohort entropy ($W=111.5$, $p<.001$), phoneme surprisal ($W=100$,
$p<.001$), word frequency ($W=115.5$, $p<.001$), word surprisal ($W=117.5$,
$p<.001$).

Figure 2: Match-mismatch accuracy values grouped for all features. A.
Accuracy values for all envelope-based features. B.  Accuracy values for
segmentation features (phoneme and word onsets) and linguistic features,
grouped per level (phoneme/word level). Group comparisons were performed using
Wilcoxon rank-sum tests. ***: $p<.001$, **: $p<.01.$

Next, we used individual participants’ match-mismatch accuracy values for all
11 features, together with age, as input to an SVM to classify the aphasia
status of participants using leave-one-out cross-validation. The SVM
classified participants belonging to either group with an accuracy of 85.42%,
an F1-score of 85.42% and an AUC of 86.19%. The SVM had a sensitivity of
92.31% and a specificity of 77.27% for aphasia. Figure 3A displays the ROC
curve. The SHAP analysis identified theta band envelope tracking as the
strongest discriminative feature (see Figure 3B).

Figure 3: SVM performance and feature importance ranking. A. ROC curve,
plotting the false positive rate as a function of the true positive rate. B.
SHAP analysis investigating the impact of each feature on the SVM’s
predictions. Features are ranked according to their impact. C. SVM performance
as a function of EEG recording length. The 9-minute mark, achieving the same
accuracy as the entire recording, is indicated in gray.

Finally, Figure 3C displays the SVM classification accuracy as a function of
the recording length. At the 9-minute mark, SVM accuracy approaches that of
the full recording and consistently surpasses the 80% accuracy and AUC
benchmark thereafter.

##  IV Discussion

We reported decreased neural tracking in IWA compared to healthy controls for
acoustic, segmentation, and linguistic speech features. The utility of neural
tracking as an effective marker for aphasia was demonstrated using an SVM
classifier which detected aphasia with high accuracy (85%) and in a time-
efficient manner (requiring 9 minutes of data).

Decreased tracking of acoustic and segmentation features in IWA compared to
healthy controls aligns with previous research on the same dataset [8, 9].
Leveraging the more complex and nonlinear architecture of CNN-based models,
our approach successfully differentiated both groups in neural tracking of
linguistic features as well, which linear models did not achieve [8]. Using
all speech features as input to an SVM, we obtained a small classification
improvement of 2% accuracy over a previous SVM classifier applied to the same
dataset, which relied on mutual information-based tracking of the speech
envelope only [9]. It is plausible that SVM performance is at ceiling level in
our sample, primarily consisting of mildly impaired IWA (27% did not score
below the clinical threshold for aphasia at time of study participation).
Alternatively, the small improvement may mean that segmentation and linguistic
features contribute only modestly over the envelope for aphasia
classification. Nevertheless, this paper introduces a substantial enhancement
compared to prior work in its clinical application potential, as CNN-based
models demonstrate generalization to unseen data and speech materials. This
eliminates the need for subject-specific training and enhances applicability
to different settings.

Neural tracking of the low-frequency envelope (delta and theta) emerged as the
most influential feature for aphasia classification (see Figure 3B). The
envelope of speech is an essential cue for speech understanding [19],
encompassing cues for processing linguistic structures of speech [10], and can
be predicted from segmentation and linguistic speech features [5]. This may
explain the relatively low contribution of higher-level segmentation and
linguistic features over the envelope for aphasia detection. Furthermore, the
present approach including segmentation and linguistic features required
slightly more data for optimal SVM performance compared to our previous
results using only the envelope as input (9 vs. 7 minutes) [9]. This suggests
that CNNs may require more input data from linguistic features, aligning with
prior findings using linear models [20].

The presented approach exhibited remarkably high sensitivity to detect aphasia
(92%), surpassing the 73% of the specific behavioral tests we administered in
our study (see Methods and [9]), suggesting a more sensitive screening tool
for aphasia. However, before implementation in clinical settings, attention
must be given to enhancing the specificity of our approach (currently at 77%)
which may require larger datasets. Furthermore, our approach must be validated
against extensive behavioral test batteries, and its ability to capture
impaired language components (i.e., specific problems with e.g. phonetic or
semantic processes) must be assessed. Additionally, testing our approach on
more severe IWA and in the acute post-stroke stage is essential, along with
validation against stroke patients without aphasia (a critique that applies to
behavioral tests as well, see [4]). Nevertheless, our work represents a
significant step towards more automatic, objective, and time-efficient
assessments of natural speech in aphasia.

##  V Conclusion

The current study presents an automated screening tool for speech processing
impairments in post-stroke aphasia using brain responses to natural speech in
a deep learning framework. By leveraging deep neural networks trained on
diverse speech materials and a large dataset of healthy individuals, which
serves as a model for intact speech processing, our approach demonstrates the
ability to generalize to unseen participants and to detect impaired speech
processing in aphasia. With its high robustness, generalizability, and time-
efficiency, our approach holds significant promise for clinical applications
and opens new perspectives to studying various other language disorders with
pre-trained CNNs.

## Acknowledgment

The authors would like to thank all participants. Thanks to Dr. Klara
Schevenels for her assistance in the recruitment process and all individuals
who helped with the data collection. Research was financially supported by the
Research Foundation Flanders (PhD grant PDC: 1S40122N, PhD grant CP: 1S49823N,
postdoc grant JV: 1290821, grant MV and TF: G0D8520N), the Luxembourg National
Research Fund (FNR, AFR-PhD project JK: 13513810) and the KU Leuven Special
Research Fund C24/18/099 (C2 project HVH and TF: 637424).

## References

  * [1] B. Gialanella, M. Bertolinelli, M. Lissi, and P. Prometti, “Predicting outcome after stroke: The role of aphasia,” _Aphasiology_ , vol. 33, no. 2, p. 122–129, 2010.
  * [2] B. N. Pasley and R. T. Knight, “Decoding speech for understanding and treating aphasia,” _Progress in Brain Research_ , vol. 207, pp. 435–456, 2013.
  * [3] H. El Hachioui, E. G. Visch-Brink, H. F. Lingsma, M. W. Van De Sandt-Koenderman, D. W. Dippel, P. J. Koudstaal, and H. A. Middelkoop, “Nonlinguistic cognitive impairment in poststroke aphasia: A prospective study,” _Neurorehabilitation and Neural Repair_ , vol. 28, no. 3, pp. 273–281, 2014.
  * [4] A. Rohde, L. Worrall, E. Godecke, R. O’Halloran, A. Farrell, and M. Massey, “Diagnosis of aphasia in stroke populations: A systematic review of language tests,” _PLoS ONE_ , vol. 13, no. 3, p. e0194143, 2018.
  * [5] M. Gillis, J. Van Canneyt, T. Francart, and J. Vanthornhout, “Neural tracking as a diagnostic tool to assess the auditory pathway,” _Hearing Research_ , p. 108607, 2022.
  * [6] M. J. Crosse, G. M. Di Liberto, A. Bednar, and E. C. Lalor, “The multivariate temporal response function (mTRF) toolbox: A MATLAB toolbox for relating neural signals to continuous stimuli,” _Frontiers in Human Neuroscience_ , vol. 10, no. NOV2016, pp. 1–14, 2016.
  * [7] A. de Cheveigné, M. Slaney, S. A. Fuglsang, and J. Hjortkjaer, “Auditory stimulus-response modeling with a match-mismatch task,” _Journal of Neural Engineering_ , vol. 18, no. 4, p. 046040, aug 2021. [Online]. Available: https://iopscience.iop.org/article/10.1088/1741-2552/abf771
  * [8] J. Kries, P. De Clercq, M. Gillis, J. Vanthornhout, R. Lemmens, T. Francart, and M. Vandermosten, “Exploring neural tracking of acoustic and linguistic speech representations in individuals with post-stroke aphasia,” _bioRxiv_ , p. 2023.03.01.530707, 2023.
  * [9] P. De Clercq, J. Kries, R. Mehraram, J. Vanthornhout, T. Francart, and M. Vandermosten, “Detecting post-stroke aphasia using eeg-based neural envelope tracking of natural speech,” _medRxiv_ , p. 2023.03.14.23287194, 2023.
  * [10] A. Giraud and D. Poeppel, “Cortical oscillations and speech processing: emerging computational principles and operations,” _Nat. Neurosci._ , vol. 15, pp. 511–517, 2012.
  * [11] C. Puffay, J. Vanthornhout, M. Gillis, B. Accou, H. V. hamme, and T. Francart, “Robust neural tracking of linguistic speech representations using a convolutional neural network,” _Journal of Neural Engineering_ , vol. 20, no. 4, p. 046040, aug 2023.
  * [12] C. Puffay, B. Accou, L. Bollens, M. J. Monesi, J. Vanthornhout, H. V. hamme, and T. Francart, “Relating eeg to continuous speech using deep neural networks: a review,” _Journal of Neural Engineering_ , vol. 20, no. 4, p. 041003, aug 2023.
  * [13] B. Accou, M. Jalilpour-Monesi, J. Montoya-Martínez, H. V. hamme, and T. Francart, “Modeling the relationship between acoustic stimulus and eeg with a dilated convolutional neural network,” _2020 28th European Signal Processing Conference (EUSIPCO)_ , pp. 1175–1179, 2021.
  * [14] C. Puffay, J. Van Canneyt, J. Vanthornhout, H. Van hamme, and T. Francart, “Relating the fundamental frequency of speech with EEG using a dilated convolutional network,” _23rd annual conference of the International Speech Communication Association (ISCA) - Interspeech 2022_ , pp. 4038–4042, 2022.
  * [15] B. Accou, M. J. Monesi, H. V. hamme, and T. Francart, “Predicting speech intelligibility from eeg in a non-linear classification paradigm*,” _Journal of Neural Engineering_ , vol. 18, no. 6, p. 066008, nov 2021. [Online]. Available: https://dx.doi.org/10.1088/1741-2552/ac33e9
  * [16] Y. Benjamini and Y. Hochberg, “Controlling the false discovery rate: a practical and powerful approach to multiple testing,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 57, no. 1, pp. 289–300, 1995.
  * [17] Gillis, Kries, M. Vandermosten, and T. Francart, “Neural tracking of linguistic and acoustic speech representations decreases with advancing age,” _NeuroImage_ , vol. 267, no. 119841, pp. 1–16, 2 2023.
  * [18] S. M. Lundberg and S.-I. Lee, “A unified approach to interpreting model predictions,” in _Advances in Neural Information Processing Systems_ , 2017.
  * [19] R. V. Shannon, F. G. Zeng, K. V, J. Wygonski, and E. M. S, “Speech recognition with primarily temporal cues,” _Science_ , vol. 270, no. 5234, pp. 303–304, 1995.
  * [20] J. Mesik and M. Wojtczak, “The effects of data quantity on performance of temporal response function analyses of natural speech processing,” _bioRxiv_ , 2022.