$a$$a$affiliationtext: Department of Statistics, University of
Washington$b$$b$affiliationtext: Seminar for Statistics, ETH
Zürichccaffiliationtext: Departments of Computing and Mathematical Sciences
and of Electrical Engineering, California Institute of Technology

# Model Selection over Partially Ordered Sets

Armeen Taeb  Peter Bühlmann  Venkat Chandrasekaran

###### Abstract

In problems such as variable selection and graph estimation, models are
characterized by Boolean logical structure such as presence or absence of a
variable or an edge. Consequently, false positive error or false negative
error can be specified as the number of variables/edges that are incorrectly
included or excluded in an estimated model. However, there are several other
problems such as ranking, clustering, and causal inference in which the
associated model classes do not admit transparent notions of false positive
and false negative errors due to the lack of an underlying Boolean logical
structure. In this paper, we present a generic approach to endow a collection
of models with partial order structure, which leads to a hierarchical
organization of model classes as well as natural analogs of false positive and
false negative errors. We describe model selection procedures that provide
false positive error control in our general setting and we illustrate their
utility with numerical experiments.

Keywords: combinatorics $|$ greedy algorithms $|$ multiple testing $|$
stability

##  1 Introduction

In data-driven approaches to scientific discovery, one is commonly faced with
the problem of model selection. Popular examples include variable selection
(which covariates predict a response?) and graph estimation (which pairs of
variables have nonzero correlation or partial correlation?). As exemplified by
these two problems, a common feature of most model selection problems in the
literature is that the collection of models is organized according to some
type of Boolean logical structure, such as the presence versus absence of a
variable or an edge. A consequence of such structure is that model complexity
can be conveniently specified as the number of attributes (variables or edges)
in a model, while false positive error or false negative error corresponds to
the number of attributes that are incorrectly included or excluded in the
model.

In many contemporary applications, models represent a far richer range of
phenomena that are not conveniently characterized via Boolean logical
structure. As a first example, suppose we are given observations of covariate-
response pairs and we wish to order the covariates based on how well they
predict a response; the collection of models is given by the set of rankings
of the covariates. Second, consider a clustering problem in which we are given
observations of a collection of variables and the goal is to group them
according to some measure of affinity, with the number of groups and the
number of variables assigned to each group not known a priori; here the model
class is given by the collection of all possible partitions of the set of
variables. Third, suppose we wish to identify causal relations underlying a
collection of variables; the model class is the set of completed partially
directed acyclic graphs. Finally, consider the blind source separation problem
in which we are given a signal expressed as an additive combination of source
signals and our objective is to identify the constituent sources, without
prior information about the number of sources or their content; here the model
class is the collection of all possible linearly independent subsets of
vectors.

In these preceding examples, we lack a systematic definition of model
complexity, false positive error, and false negative error due to the absence
of Boolean logical structure in each collection of models. In particular, in
the first three examples, valid models are characterized by structural
properties such as transitivity, set partitioning, and graph acyclicity,
respectively; these properties are global in nature and are not concisely
modeled via separable and local characteristics such as an attribute (a
variable or edge) being included in a model independently of other attributes.
In the fourth example of blind source separation, false positive and false
negative errors should not be defined merely via the inclusion or exclusion of
true source vectors in an estimated set but should instead consider the degree
of alignment between the estimated and true sources, which again speaks to the
lack of a natural Boolean logical structure underlying the associated model
class.

As a concrete illustration of the inappropriateness of Boolean logical
structure for clustering, consider three items $a,b,c$, with the null model
given by the three clusters $\\{a\\},\\{b\\},\\{c\\}$, the true model by the
two clusters $\\{a,b\\},\\{c\\}$, and the estimated model by the single
cluster $\\{a,b,c\\}$. An incorrect perspective grounded in Boolean logical
structure suggests a false positive error of two, with the ‘false discoveries’
being that $c$ is in the same cluster as $a$ and as $b$. However, accounting
for set partition structure yields the more accurate false positive error
value of one as $a$ and $b$ are in the same cluster in the true and estimated
models; hence, including $c$ in the same cluster as $\\{a,b\\}$ should only
incur one false discovery.

While the preceding four problems have been studied extensively, the
associated methods do not systematically control false positive error as this
quantity is not formally defined. Selection procedures that yield models with
small false positive error play an important role in data-driven methods for
gathering evidence, rooted in the empirical philosophy and statistical testing
foundations of falsification of theories and hypotheses [17, 7, 18].

###  1.1 Our Contributions

We begin in Section 2 by describing how collections of models may be endowed
with the structure of a partially ordered set (poset). Posets are relations
that satisfy reflexivity, transitivity, and antisymmetry, and they facilitate
a hierarchical organization of a set of models that leads to a natural
definition of model complexity. Building on this framework, we develop an
axiomatic approach to defining functions over poset element pairs for
evaluating similarity. This yields generalizations of well-known measures such
as family-wise error and false discovery rate to an array of model selection
problems in the context of ranking, causal inference, multiple change-point
estimation, clustering, multi-sample testing, and blind source separation. In
Section 3, we describe two generic model selection procedures that search over
poset elements in a greedy fashion and that provide false discovery control in
discrete model posets. The first method is based on subsampling and model
averaging and it builds on the idea of stability selection [14, 23] for the
variable selection problem, while the second method considers a sequence of
hypothesis tests between models of growing complexity. With both these
methods, the combinatorial properties of a model poset play a prominent role
in determining computational and statistical efficiency. Proofs of the
theorems of Section 3 are provided in Section 5. In Section 4 we provide
numerical illustration via experiments on synthetic and real data. The code
for implementing our methods is available at
https://github.com/armeentaeb/model-selection-over-posets.

###  1.2 Related Work

Classic approaches to model selection such as the AIC and BIC assess and
penalize model complexity by counting the number of attributes included in a
model [1, 22]. More generally, such complexity measures facilitate a
hierarchical organization of model classes, and this perspective is prevalent
throughout much of the model selection literature [9, 13, 29, 19, 8, 2].
However, these complexity measures rely on a Boolean logical structure
underlying a collection of models, and are therefore not well-suited to model
classes that are not characterized in this manner. The poset formalism
presented in this paper is sufficiently flexible to facilitate model selection
over model classes that are more complex than those characterized by Boolean
logical structure (such as the illustration presented previously with
clustering, see also Example 2), while being sufficiently structured to permit
precise definitions of model complexity as well as false positive and false
negative errors.

##  2 Poset Framework for Model Selection

We begin by describing how collections of models arising in various
applications may be organized as posets. Next, we present approaches to endow
poset-structured models with suitable notions of true and false discoveries.

###  2.1 Model Classes as Posets

We begin with some basics of posets. A _poset_ $(\mathcal{L},\preceq)$ is a
collection $\mathcal{L}$ of elements and a relation $\preceq$ that is
reflexive ($x\preceq x,~{}\forall x\in\mathcal{L}$), transitive ($x\preceq
y,y\preceq z\Rightarrow x\preceq z,~{}\forall x,y,z\in\mathcal{L}$), and anti-
symmetric ($x\preceq y,y\preceq x\Rightarrow x=y,~{}\forall
x,y\in\mathcal{L}$). An element $y\in\mathcal{L}$ _covers_ $x\in\mathcal{L}$
if $x\preceq y$, $x\neq y$, and there is no
$z\in\mathcal{L}\backslash\\{x,y\\}$ with $x\preceq z\preceq y$; we call such
$(x,y)$ a _covering pair_. A _path_ from $x_{1}\in\mathcal{L}$ to
$x_{k}\in\mathcal{L}$ is a sequence $(x_{1},\dots,x_{k})$ with
$x_{2},\dots,x_{k-1}\in\mathcal{L}$ such that $x_{i}$ covers $x_{i-1}$ for
each $i=2,\dots,k$. Throughout this paper, we focus on posets in which there
is a _least element_ , i.e., an element $x_{\text{least}}\in\mathcal{L}$ such
that $x_{\text{least}}\preceq y$ for all $y\in\mathcal{L}$; such least
elements are necessarily unique. Finally, a poset is _graded_ if there exists
a function $\mathrm{rank}(\cdot)$ mapping poset elements to the nonnegative
integers such that the rank of the least element is $0$ and
$\mathrm{rank}(y)=\mathrm{rank}(x)+1$ for $y\in\mathcal{L}$ that covers
$x\in\mathcal{L}$. In graded posets with least elements, each path from the
least element to any $x\in\mathcal{L}$ has length equal to $\mathrm{rank}(x)$.
Posets are depicted visually using Hasse diagrams in which a directed arrow is
drawn from $x\in\mathcal{L}$ to any $y\in\mathcal{L}$ that covers $x$.

Posets offer an excellent framework to formulate model selection problems as
model classes in many applications possess rich partial order structures. In
particular, the poset-theoretic quantities in the preceding paragraph have
natural counterparts in the context of model selection – the least element
corresponds to the ‘null’ model that represents no discoveries, the relation
$\preceq$ specifies a notion of containment between simpler and more complex
models, and the rank function serves as a measure of model complexity that
respects the underlying containment relation. We present several concrete
illustrations next; Figure 1 presents Hasse diagrams associated with several
of these examples.

(a) Variable selection (b) Clustering (c) Multisample testing (d) Causal
inference (e) Partial ranking (f) Total ranking

Figure 1: Hasse diagrams for a) variable selection with $3$ variables (Example
1); b) clustering $4$ variables (Example 2); c) multisample testing with $4$
samples (Example 3); d) causal inference with $3$ variables (Example 4); e)
partial ranking of $3$ items (Example 6); and f) total ranking of $3$ items
(Example 7).

######  Example 1 (Variable selection).

As a warm-up, consider the variable selection problem of selecting which of
$p$ variables influence a response. The poset here is the collection of all
subsets of $\\{1,\dots,p\\}$ ordered by set inclusion, the least element is
given by the empty set, and the rank of a subset is its cardinality. This
poset is called the _Boolean poset_ [24].

######  Example 2 (Clustering).

Suppose we wish to group a collection of $p$ variables based on a given notion
of similarity. The poset here is the collection of all partitions of
$\\{1,\dots,p\\}$ ordered by refinement, the least element is given by $p$
groups each consisting of one variable, and the rank of a partition is equal
to $p$ minus the number of groups. Thus, higher-rank elements correspond to
models specified by a small number of clusters. This poset is called the
_partition poset_ [24].

######  Example 3 (Multisample testing).

As a generalization of the classic two-sample testing problem, consider the
task of grouping $p$ samples with the objective that samples in a group come
from the same distribution. Although this problem is closely related to the
preceding clustering problem, it is more natural for the underlying poset here
to be the reverse of the partition poset that is formed by reversing the order
relation of the partition poset, i.e., the poset is the collection of all
partitions of $\\{1,\dots,p\\}$ ordered by coarsening. With this reverse
ordering, the least element corresponds to all $p$ samples belonging to the
same group (i.e., coming from the same distribution), which generalizes the
usual null hypothesis in two-sample testing. The rank of a partition is equal
to the number of groups minus one. Thus, higher-rank elements correspond to
the $p$ samples arising from many distinct distributions.

######  Example 4 (Causal structure learning).

Causal associations among a collection of variables are often characterized by
a directed acyclic graph (DAG), namely a graph with directed edges and no
(directed) cycles, in which the nodes index the variables. Causal structure
learning entails inferring this DAG from observations of the variables. The
structure of a DAG specifies a causal model via conditional independence
relations among the variables, with denser DAGs encoding fewer conditional
independencies in comparison with sparser DAGs. (See [5] for details on how
the structure of a DAG encodes conditional independence relations; here we
describe only those aspects that pertain to a poset formulation to organize
the collection of all causal models based on graph structure.) Distinct DAGs
can specify the same set of conditional independence relations, and these are
called Markov equivalent DAGs. We introduce some terminology to characterize
Markov equivalent DAGs. The _skeleton_ of a DAG is the undirected graph
obtained by making all the edges undirected. A _$v$ -structure_ is a set of
three nodes $x,y,z$ such that there are directed edges from $x$ to $z$ and
from $y$ to $z$, and there is no edge between $x$ and $y$. Two DAGs are Markov
equivalent if and only if they have the same skeleton and the same collection
of $v$-structures [26]. A Markov equivalence class of DAGs can be described by
a completed partially directed acyclic graph (CPDAG), which is a graph
consisting of both directed and undirected edges. A CPDAG has a directed edge
from a node $x$ to a node $y$ if and only if this directed edge is present in
every DAG in the associated Markov equivalence class. A CPDAG has an
undirected edge between nodes $x$ and $y$ if the corresponding Markov
equivalence class contains a DAG with a directed edge from $x$ to $y$ and a
DAG with a directed edge from $y$ to $x$. One can check that the total number
of edges in a CPDAG (directed plus undirected) is equal to the number of edges
in any DAG in the associated Markov equivalence class. The collection of
CPDAGs on $p$ variables may be viewed as a poset ordered by inclusion – CPDAGs
$\mathcal{C}^{(1)},\mathcal{C}^{(2)}$ satisfy
$\mathcal{C}^{(1)}\preceq\mathcal{C}^{(2)}$ if and only if there exist DAGs
$\mathcal{G}^{(1)},\mathcal{G}^{(2)}$ in the respective Markov equivalence
classes such that $\mathcal{G}^{(1)}$ is a directed subgraph of
$\mathcal{G}^{(2)}$. In other words,
$\mathcal{C}^{(1)}\preceq\mathcal{C}^{(2)}$ if and only if all the conditional
independencies encoded by $\mathcal{C}^{(2)}$ are also encoded by
$\mathcal{C}^{(1)}$, or equivalently that all the conditional dependencies
encoded by $\mathcal{C}^{(1)}$ are also encoded by $\mathcal{C}^{(2)}$. The
least element is given by the CPDAG with no edges, and the rank function is
equal to the number of edges. Higher-rank elements in this poset correspond to
causal models exhibiting more conditional dependence relations. In some causal
inference contexts, it may be more natural to view the fully connected CPDAG
as the null model. Our framework can accommodate this perspective by reversing
the preceding model poset. Specifically, in this reversed poset, the partial
order is given by inclusion of conditional independencies, the least element
is the fully connected CPDAG, and the rank function is $p(p-1)/2$ minus the
number of edges. Higher-rank elements in this poset correspond to causal
models exhibiting more conditional independence relations.

######  Example 5 (Multiple changepoint estimation).

Consider the problem of detecting changepoints in a multivariate time series.
Specifically, we observe $p$ signals each for time instances $t=0,\dots,T-1$,
each signal consists of at most one change (e.g., a change in the distribution
or dynamics underlying the signal observations), and the objective is to
identify these changes. We denote changepoints via vectors
$x=(x_{1},\dots,x_{p})\in\\{0,\dots,T\\}^{p}$, with $x_{i}$ denoting the time
index when a change occurs in the $i$’th signal and $x_{i}=T$ corresponding to
no change occurring. The poset here is the set $\\{0,1,\dots,T\\}^{p}$ ordered
such that $x\preceq y$ if and only if $x_{j}\geq y_{j}$ for all $j=1,\dots,p$,
the least element is $(T,\dots,T)$, and the rank of an element is $p\cdot T$
minus the sum of the coordinates. Higher-rank elements correspond to
changepoint estimates in which the changes occur early. This poset is the
reverse of the (bounded) _integer poset_ [24] with the product order.

######  Example 6 (Partial ranking).

We seek a ranking of a finite set of items given noisy observations (e.g.,
pairwise comparisons), and we allow some pairs of items to be declared as
incomparable. Such a _partial ranking_ of the elements of a finite set $S$
corresponds to a _strict partial order_ on $S$, i.e., a relation $\mathcal{R}$
that is irreflexive $((a,a)\notin\mathcal{R},~{}\forall a\in S)$, asymmetric
$((a,b)\in\mathcal{R}\Rightarrow(b,a)\notin\mathcal{R},~{}\forall a,b\in S)$,
and transitive; if an element of $S$ does not appear in $\mathcal{R}$, then
that element is incomparable to any of the other elements of $S$ in the
associated partial ranking. The poset here is the collection of strict partial
orders on $S$ ordered by inclusion, the least element is the empty set, and
the rank of a partial ranking is the cardinality of the associated relation.
Thus, higher-rank elements correspond to partial rankings that compare many of
the covariates.

######  Example 7 (Total ranking).

We again wish to rank a finite collection of items but now we seek a total
ranking that provides an ordered list of all the items. The setting is that we
are given a total ranking that represents our current state of knowledge
(i.e., a ‘null model’) as well as a new set of noisy observations, and the
goal is to identify a total ranking that represents an update of the null
model to reflect the new information. Each _total ranking_ of the elements of
a finite set $S$ corresponds to a one-to-one function from $S$ to the integers
$\\{1,\dots,|S|\\}$. Let $\pi_{\mathrm{null}}$ be the function that describes
the null ranking. A convenient way to compare total rankings and to define a
poset structure over them is via the notion of an _inversion set_. For any
total ranking specified by a function $\pi$, the associated inversion set
(with respect to the null ranking $\pi_{\mathrm{null}}$) is defined as
$\mathrm{inv}(\pi;\pi_{\mathrm{null}}):=\\{(x,y)\in S\times
S~{}|~{}\pi_{\mathrm{null}}(x)<\pi_{\mathrm{null}}(y),~{}\pi(x)>\pi(y)\\}$.
The poset here (with respect to a given null ranking $\pi_{\mathrm{null}}$) is
the collection of total rankings on $S$ ordered by inclusion of the associated
inversion sets, the least element is the null ranking $\pi_{\mathrm{null}}$,
and the rank of a total ranking is the cardinality of the associated inversion
set; this rank function is also equal to the Kendall tau distance between a
total ranking and $\pi_{\mathrm{null}}$. Thus, higher-rank elements are given
by total rankings that depart significantly from the null ranking
$\pi_{\mathrm{null}}$. This poset is called the _permutation poset_ [24].

######  Example 8 (Subspace estimation).

The task is to estimate a subspace in $\mathbb{R}^{p}$ given noisy
observations of points in the subspace. The poset is the collection of
subspaces in $\mathbb{R}^{p}$ ordered by inclusion, the least element is the
subspace $\\{0\\}$, and the rank of a subspace is its dimension. This poset is
called the _subspace poset_.

######  Example 9 (Blind source separation).

We are given a signal in $\mathbb{R}^{p}$ that is expressed as a linear
combination of some unknown source signals and the goal is to estimate these
sources. The poset here is the collection of linearly independent subsets of
unit-norm vectors in $\mathbb{R}^{p}$ ordered by inclusion, the least element
is the empty set, and the rank of a linearly independent subset is equal to
the cardinality of the subset.

With respect to formalizing the notion of false positive and false negative
errors, Example 1 is prominently considered in the literature, while Examples
3 and 5 are multivariate generalizations of previously studied cases [10, 12].
Finally, Example 8 was studied in [25], although that treatment proceeded from
a geometric perspective rather than the order-theoretic approach presented in
this paper. With the exception of Example 1, none of the other examples permit
a natural formulation within the traditional multiple testing paradigm due to
the lack of a Boolean logical structure underlying the associated model
classes. Moreover, Examples 8-9 are model classes consisting of infinitely
many elements. Nonetheless, we describe in the sequel how the poset formalism
enables a systematic and unified framework for formulating model selection in
all of the examples above.

###  2.2 Evaluating True and False Discoveries

To assess the extent to which an estimated model signifies discoveries about
the true model, we describe next a general approach to quantify the similarity
between poset elements in a manner that respects partial order structure.

######  Definition 1 (similarity valuation).

Let $(\mathcal{L},\preceq,\mathrm{rank}(\cdot))$ be a graded poset. A function
$\rho:\mathcal{L}\times\mathcal{L}\to\mathbb{R}$ that is symmetric, i.e.,
$\rho(x,y)=\rho(y,x)$ for all $x,y\in\mathcal{L}$, is called a _similarity
valuation_ over $\mathcal{L}$ if:

  * •

$0\leq\rho(x,y)\leq\min\\{\mathrm{rank}(x),\mathrm{rank}(y)\\}$ for all
$x,y\in\mathcal{L}$,

  * •

$\rho(x,y)\leq\rho(z,y)$ for all $x\preceq z$,

  * •

$\rho(x,y)=\mathrm{rank}(x)$ if and only if $x\preceq y$.

######  Remark 1.

The term ‘valuation’ is often used in the order-theory literature [24] to
denote functions on posets that respect the underlying partial order
structure, and we use it in our context for the same reason.

In the sequel, we describe similarity valuations for the various model posets
discussed previously. The conditions above make similarity valuations well-
suited for quantifying the amount of discovery in an estimated model with
respect to a true model. The first condition states that the amount of
discovery must be bounded above by the complexities of the true and estimated
models (which are specified by the rank function). The second condition
requires similarity valuations to respect partial order structure so that more
complex models do not yield less discovery than less complex ones. The final
condition expresses the desirable property that the amount of discovery
contained in an estimated model is equal to the complexity of that model if
and only if it is ‘contained in’ the true model. With these properties, we
obtain the following analogs of true and false discoveries and of related
quantities such as false discovery proportion.

######  Definition 2 (true and false discoveries).

Let $(\mathcal{L},\preceq,\mathrm{rank}(\cdot))$ be a graded poset and let
$\rho$ be a similarity valuation on $\mathcal{L}$. Letting
$x^{\star}\in\mathcal{L}$ be a true model and $\hat{x}\in\mathcal{L}$ be an
estimate, the _true discovery_ , the _false discovery_ , and the _false
discovery proportion_ are, respectively, defined as follows:

$\displaystyle\mathrm{TD}(\hat{x},x^{\star}):=\rho(\hat{x},x^{\star}),$
$\displaystyle\mathrm{FD}(\hat{x},x^{\star}):=\mathrm{rank}(\hat{x})-\rho(\hat{x},x^{\star})=\mathrm{rank}(\hat{x})-\mathrm{TD}(\hat{x},x^{\star}),$
$\displaystyle\mathrm{FDP}(\hat{x},x^{\star}):=\frac{\mathrm{rank}(\hat{x})-\rho(\hat{x},x^{\star})}{\mathrm{rank}(\hat{x})}=\frac{\mathrm{FD}(\hat{x},x^{\star})}{\mathrm{rank}(\hat{x})}.$

With these definitions, we articulate our model selection objective more
precisely:

Goal: _identify the largest rank model subject to control in expectation or in
probability on false discovery (proportion)._

This objective is akin to seeking the largest amount of discovery subject to
control on false discovery (rate). The data available to carry out model
selection vary across our examples; in Section 3 we describe methods to obtain
false discovery control guarantees in various settings.

To carry out this program, a central question is the choice of a suitable
similarity valuation for a graded model poset. Indeed, it is unclear whether
there always exists a similarity valuation for any graded model poset
$(\mathcal{L},\preceq,\mathrm{rank}(\cdot))$. To address this question,
consider the following function for $x,y\in\mathcal{L}$:

$\rho_{\text{meet}}(x,y):=\max_{z\preceq x,z\preceq y}~{}\mathrm{rank}(z).$
(1)

######  Remark 2.

In order theory, a poset $(\mathcal{L},\preceq)$ is said to possess a _meet_
if for each $x,y\in\mathcal{L}$ there exists a $z\in\mathcal{L}$ satisfying
$(i)$ $z\preceq x,z\preceq y$ and $(ii)$ for any $w\in\mathcal{L}$ with
$w\preceq x,w\preceq y$, we have $w\preceq z$; such a $z$ is called the meet
of $x,y$ and posets that possess a meet are called _meet semi-lattices_.
Except for the poset in Example 4 on causal structure learning, the posets in
the other examples are meet semi-lattices (see SI Appendix, Section 1). The
subscript ‘meet’ in (1) signifies that $\rho_{\text{meet}}$ is the rank of the
meet for meet semi-lattices, although $\rho_{\text{meet}}$ is well-defined
even if $(\mathcal{L},\preceq)$ is not a meet semi-lattice.

One can check that $\rho_{\text{meet}}$ is a similarity valuation on any
graded poset $(\mathcal{L},\preceq,\mathrm{rank}(\cdot))$; see SI Appendix,
Section 2 for a proof. For Example 1 on variable selection,
$\rho_{\text{meet}}$ has the desirable property that it reduces to the number
of common variables in two models; thus, the general model selection goal
formulated above reduces to the usual problem of maximizing the number of
selected variables subject to control on the number of selected variables that
are null. Next, we describe the model selection problems we obtain in Examples
2-6 with $\rho_{\text{meet}}$ as the choice of similarity valuation.

In Example 2 on clustering, the value of $\rho_{\text{meet}}$ for two
partitions of $p$ variables is equal to $p$ minus the number of groups in the
coarsest common refinement of the partitions. The model selection problem is
that of partitioning the variables into the smallest number of groups subject
to control on the additional number of groups in the coarsest common
refinement of the estimated and true partitions compared to the number of
groups in the estimated partition.

Recall that the poset in Example 3 on multisample testing is the reverse of
the poset in Example 2; thus, many of the notions from the preceding paragraph
are appropriately ‘reversed’ in Example 3. In particular, the value of
$\rho_{\text{meet}}$ in Example 3 for two partitions of $p$ samples is equal
to the number of groups in the finest common coarsening of the partitions. The
model selection problem entails partitioning the samples into the largest
number of groups subject to control on the additional number of groups in the
estimated partition compared to the number of groups in the finest common
coarsening of the estimated and true partitions.

In Example 4 on causal structure learning, the value of $\rho_{\text{meet}}$
for two CPDAGs $\mathcal{C}^{(1)},\mathcal{C}^{(2)}$ is equal to the maximum
number of edges in a CPDAG that encodes all the conditional independencies of
$\mathcal{C}^{(1)}$ and of $\mathcal{C}^{(2)}$. The model selection task is
then to identify the CPDAG with the largest number of edges subject to control
on the additional number of edges in the estimated CPDAG compared to the
densest CPDAG that encodes all the conditional independence relationships in
both the true and estimated CPDAGs.

In Example 5 on multiple changepoint estimation, suppose
$x,y\in\\{0,\dots,T\\}^{p}$ are vectors of time indices specifying
changepoints in $p$ signals. We have that
$\rho_{\text{meet}}(x,y)\allowbreak=p\cdot
T-\allowbreak\sum_{i=1}^{p}\max\\{x_{i},y_{i}\\}$. The model selection problem
entails identifying changes as quickly as possible subject to control on early
detection of changes (i.e., declaring changes before they occur); this is a
multivariate generalization of the classic quickest change detection problem
[12].

In Example 6 on partial ranking, the value of $\rho_{\text{meet}}$ for two
partial rankings is equal to the cardinality of the intersection of the
associated relations, i.e., the number of common comparisons in the two
partial rankings. The associated model selection problem is that of
identifying a partial ranking with the largest number of comparisons (i.e.,
the associated relation must have large cardinality) subject to control on the
number of comparisons in the estimated partial ranking that are not in the
true partial ranking.

In Examples 1-6, the function $\rho_{\text{meet}}$ of (1) provides a
convenient way to assess the amount of discovery in an estimated model with
respect to a true model, thereby yielding natural formulations for model
selection. However, in Examples 7-9, $\rho_{\text{meet}}$ has some undesirable
features.

Consider first the setup in Example 7 on total ranking for the set
$S=\\{a,b,c\\}$ with the null model given by the ranking
$\pi_{\text{null}}(a)=1,\pi_{\text{null}}(b)=2,\pi_{\text{null}}(b)=3$, the
true model given by the ranking
$\pi^{\star}(a)=3,\pi^{\star}(b)=1,\pi^{\star}(c)=2$ (Hasse diagram shown in
Figure 1), and the estimated ranking given by
$\hat{\pi}(a)=2,\hat{\pi}(b)=3,\hat{\pi}(c)=1$. In this case, one can see from
Figure 1 that $\rho_{\text{meet}}(\hat{\pi},\pi^{\star})=0$, which suggests
that no discovery is made. On the other hand, the inversion sets of these
rankings are given by
$\mathrm{inv}(\pi^{\star};\pi_{\text{null}})=\\{(a,b),(a,c)\\}$ and
$\mathrm{inv}(\hat{\pi};\pi_{\text{null}})=\\{(a,c),(b,c)\\}$, and the element
$(a,c)$ is common to both inversion sets as the fact that item $c$ is ranked
higher than item $a$ in the true model has been discovered in the estimated
model; this reasoning suggests that a positive quantity would be a more
appropriate value for the similarity valuation between $\hat{\pi}$ and
$\pi^{\star}$. The key issue is that
$\mathrm{inv}(\pi^{\star};\pi_{\text{null}})\cap\mathrm{inv}(\hat{\pi};\pi_{\text{null}})$
is not an inversion set of any total ranking, but this intersection still
carries valuable information about true discoveries made in $\hat{\pi}$ about
$\pi^{\star}$. However, the similarity valuation $\rho_{\text{meet}}$ only
considers subsets of
$\mathrm{inv}(\pi^{\star};\pi_{\text{null}})\cap\mathrm{inv}(\hat{\pi};\pi_{\text{null}})$
that correspond to inversion sets of total rankings as the maximization in (1)
is constrained to be over poset elements. Motivated by this discussion, we
employ the following similarity valuation in Example 7 for total rankings
$\pi,\tilde{\pi}$ (with respect to a null model $\pi_{\text{null}}$):

$\rho_{\text{total-
ranking}}(\pi,\tilde{\pi})=|\mathrm{inv}(\pi;\pi_{\text{null}})\cap\mathrm{inv}(\tilde{\pi};\pi_{\text{null}})|.$
(2)

With this similarity valuation, the model selection problem reduces to
identifying a total ranking with the largest inversion set (with respect to
$\pi_{\text{null}}$) subject to control on the number of comparisons in the
inversion set of the estimated total ranking that are not in the inversion set
of the true total ranking.

Next, in Example 8, $\rho_{\text{meet}}(\hat{x},x^{\star})$ is equal to the
dimension of the intersection of the subspaces $\hat{x},x^{\star}$. When these
subspaces have small dimensions, for example, $\rho_{\text{meet}}$ generically
equals zero regardless of the angle between the subspaces; in words,
$\rho_{\text{meet}}$ does not consider the smooth structure underlying the
collection of subspaces. As discussed in [25], a more suitable measure of
similarity is the sum of the squares of the cosines of the principal angles
between the subspaces, which is expressed as follows using projection matrices
onto subspaces $\mathcal{U},\tilde{\mathcal{U}}$:

$\rho_{\text{subspace}}(\mathcal{U},\tilde{\mathcal{U}})=\mathrm{trace}(\mathcal{P}_{\mathcal{U}}\mathcal{P}_{\tilde{\mathcal{U}}}).$
(3)

The model selection task is to identify the largest-dimensional subspace
subject to control on the sum of the squares of the cosines of the principal
angles between the estimated subspace and the orthogonal complement of the
true subspace.

Finally, $\rho_{\text{meet}}$ is inadequate as a similarity valuation in
Example 9 for the same reasons as in Example 8 due to the underlying smooth
structure, and we propose here a more appropriate alternative. Given
$B\in\mathbb{R}^{p\times k},~{}\tilde{B}\in\mathbb{R}^{p\times\ell}$ (these
matrices have unit-norm and linearly independent columns representing source
signals), suppose without loss of generality that $k\leq\ell$ (due to the
symmetry of similarity valuations) and let $\text{Perm}(\ell)$ be the
collection of bijections on $\\{1,\dots,\ell\\}$. With this notation, consider
the following similarity valuation:

$\rho_{\text{source-
separation}}(B,\tilde{B})=\max_{\sigma\in\text{Perm}(\ell)}\sum_{i=1}^{k}(B^{T}\tilde{B})_{i,\sigma(i)}^{2}.$
(4)

This valuation is better suited to quantify the degree of alignment between
two collections of vectors in source separation than $\rho_{\text{meet}}$.
Model selection entails identifying the largest collection of source vectors
subject to control on the difference in the number of estimated source vectors
and the alignment between the true and estimated source vectors as evaluated
by $\rho_{\text{source-separation}}$.

Table 1 summarizes our discussion of the various model posets and their
associated similarity valuations. In conclusion, while $\rho_{\text{meet}}$ is
a similarity valuation for any model poset, it is not always the most natural
choice, and identifying a suitable similarity valuation that captures the
essential features of an application is key to properly formulating a model
selection problem. This situation is not unlike the selection of an
appropriate loss function in point estimation – while there exist many
candidates that are mathematically valid, the utility of an estimation
procedure in the context of a problem domain depends critically on a well-
chosen loss.

problem domain | models |  least element (i.e. global null) |  partial order |  rank (i.e. model complexity) |  similarity valuation (i.e. true discoveries)
---|---|---|---|---|---
_ variable selection _ |  subsets of $\\{1,\dots,p\\}$ | $\emptyset$ |  inclusion of subsets |  cardinality of subset |  subsets $x,\tilde{x}$; $\rho(x,\tilde{x})=|x\cap\tilde{x}|$
_clustering_ |  partitions of $\\{1,\dots,p\\}$ | $\\{1\\},\\{2\\},\dots,\\{p\\}$ |  refinement of partition | $p-\\#\text{groups}$ |  partitions $x,\tilde{x}$; $\rho(x,\tilde{x})=p-\\#$ groups in coarsest common refinement
_ multisample testing _ |  partitions of $\\{1,\dots,p\\}$ | $\\{1,2,\dots,p\\}$ |  coarsening of partition | $\\#\text{groups}$ |  partitions $x,\tilde{x}$; $\rho(x,\tilde{x})=\\#\text{ groups in}$ finest common coarsening
_ causal structure learning _ |  completed partially directed acyclic graphs (CPDAG) on a set of variables |  CPDAG with no edges |  inclusion of conditional dependencies encoded by CPDAGs | $\\#\text{edges}$ |  $\text{CPDAGs }\mathcal{C},\tilde{\mathcal{C}}$; $\rho(\mathcal{C},\tilde{\mathcal{C}})=\\#\text{edges in }$ densest CPDAG encoding conditional independencies of both $\mathcal{C},\tilde{\mathcal{C}}$
_ multiple changepoint _ |  elements of $\\{0,\dots,T\\}^{p}$ | $(T,T,\dots,T)$ |  entrywise reverse ordering |  $p\cdot T$ minus sum of entries |  changepoint vectors $x,\tilde{x}$; $\rho(x,\tilde{x})=p\cdot T-\sum_{i}\max\\{x_{i},\tilde{x}_{i}\\}$
_ partial ranking _ |  relations specified by strict partial orders on a set of items | $\emptyset$ |  inclusion of sets specifying relations |  cardinality of set specifying relation |  sets $\mathcal{R},\tilde{\mathcal{R}}$ specifying relations; $\rho(\mathcal{R},\tilde{\mathcal{R}})=|\mathcal{R}\cap\tilde{\mathcal{R}}|$
_ total ranking _ |  total orders on a set of items |  base ranking $\pi_{\text{null}}$ |  inclusion of inversion sets w.r.t. $\pi_{\text{null}}$ |  cardinality of inversion set w.r.t. $\pi_{\text{null}}$ |  total orders $\pi,\tilde{\pi}$; $\rho(\pi,\tilde{\pi})=|\mathrm{inv}(\pi;{\pi}_{\text{null}})\cap\mathrm{inv}(\tilde{\pi};{\pi}_{\text{null}})|$
_ subspace estimation _ | subspaces in $\mathbb{R}^{p}$ | $\\{0\\}$ |  inclusion of subspaces |  dimension of subspace |  subspaces $\mathcal{U},\tilde{\mathcal{U}}$; $\rho(\mathcal{U},\tilde{\mathcal{U}})=\mathrm{trace}(\mathcal{P}_{\mathcal{U}}\mathcal{P}_{\tilde{\mathcal{U}}})$
_ blind source separation _ |  linearly independent subsets of $\mathbb{R}^{p}$ | $\emptyset$ |  inclusion of subsets |  cardinality of subset |  subsets given by columns of $B\in\mathbb{R}^{p\times k},\tilde{B}\in\mathbb{R}^{p\times\ell},k\leq\ell$; $\rho(B,\tilde{B})=\max_{\sigma\in\text{Perm}(\ell)}\sum_{i=1}^{k}(B^{T}\tilde{B})_{i,\sigma(i)}^{2}$

Table 1: Problem classes and associated characterization of model selection
via posets.

##  3 False Discovery Control over Posets

In this section, we turn our attention to the task of identifying models of
large rank that provide false discovery control. We begin in Section 3.1 with
a general greedy strategy for poset search that facilitates the design of
model selection procedures, and we specialize this framework to specific
approaches in Sections 3.2 and 3.3. Some of the discussion in Section 3.1 is
relevant for all of the posets in Examples 1-9, while the methodology
presented in Sections 3.2-3.3 is applicable to general discrete posets with
integer-valued similarity valuations such as in Examples 1-7. Along the way,
we remark on some of the challenges that arise in the two continuous cases of
Examples 8-9.

###  3.1 Greedy Approaches to Model Selection

To make progress on the problem of identifying large rank models that provide
control on false discovery, we begin by noting that the false discovery
$\mathrm{FD}(\hat{x},x^{\star})$ in an estimated model $\hat{x}$ with respect
to a true model $x^{\star}$ may be expressed as the following telescoping sum
for any path $(x_{0},x_{1},\dots,x_{k-1},x_{k})$ with $x_{0}$ being the least
element $x_{\text{least}}$ and $x_{k}=\hat{x}$:

$\mathrm{FD}(\hat{x},x^{\star})=\sum_{i=1}^{k}1-[\rho(x_{i},x^{\star})-\rho(x_{i-1},x^{\star})].$
(5)

The term $1-[\rho(x_{i},x^{\star})-\rho(x_{i-1},x^{\star})]$ may be
interpreted as the “additional false discovery” incurred by the model $x_{i}$
relative to the model $x_{i-1}$. The above decomposition of false discovery in
terms of a path from the least element to an estimated model suggests a
natural approach for model selection. In particular, we observe that a
sufficient condition for $\mathrm{FD}(\hat{x},x^{\star})$ to be small is for
each term in the above sum to be small. Thus, we will greedily grow a path
starting from the least element $x_{0}=x_{\text{least}}$ by adding one element
$x_{i}$ at a time such that each $(x_{i-1},x_{i})$ is a covering pair and each
$1-[\rho(x_{i},x^{\star})-\rho(x_{i-1},x^{\star})]$ is small. We continue this
process until we can no longer guarantee that
$1-[\rho(x_{i},x^{\star})-\rho(x_{i-1},x^{\star})]$ is small.

For such a procedure to be fruitful, we require some data-driven method to
bound $1-[\rho(x_{i},x^{\star})-\rho(x_{i-1},x^{\star})]$ as the true model
$x^{\star}$ is not known. Our objective, therefore, is to design a data-
dependent function $\Psi:\\{(a,b)~{}|~{}b\text{ covers }a\text{ in
}\mathcal{L}\\}\to[0,1]$ that takes as input covering pairs and outputs a
number in the interval $[0,1]$, and further satisfies the property that
$\Psi(u,v)$ being small is a sufficient condition for
$1-[\rho(v,x^{\star})-\rho(u,x^{\star})]$ to be small (in expectation or in
probability). Given such a function, we grow a path using the greedy strategy
outlined above by identifying at each step a covering pair that minimizes
$\Psi$. Algorithm 1 provides the details. In Sections 3.2 and 3.3, we present
two approaches for designing suitable functions $\Psi$: one based on a notion
of stability and the other based on testing. Proofs that both these methods
control for false discoveries are presented in Section 5.

Algorithm 1 Greedy sequential algorithm for model selection

1:Input: poset $\mathcal{L}$, threshold $\alpha\in[0,1]$; data-dependent
function $\Psi:\\{(a,b)~{}|~{}b\text{ covers }a\text{ in
}\mathcal{L}\\}\to[0,1]$

2:Greedy selection: Set $u=x_{\text{least}}$ and perform:

  * (a)

find $v_{\text{opt}}\in\operatorname*{argmin}_{\\{(u,v)~{}|~{}v\text{ covers
}u\text{ in }\mathcal{L}\\}}{\Psi}(u,v)$.

  * (b)

if ${\Psi}(u,v_{\text{opt}})\leq\alpha$, set $u=v_{\text{opt}}$ and repeat
steps (2a-2b). Otherwise, stop.

3:Output: return $\hat{x}=u$

In designing a suitable function $\Psi$ so that
$1-(\rho(v,x^{\star})-\rho(u,x^{\star}))$ is small (in expectation or in
probability) whenever $\Psi(u,v)$ is small, we note that the examples
presented in Section 2 exhibit an important invariance. Specifically, in each
example there are distinct covering pairs $(u,v)$ and
$(u^{\prime},v^{\prime})$ such that
$1-[\rho(v,x^{\star})-\rho(u,x^{\star})]=1-[\rho(v^{\prime},x^{\star})-\rho(u^{\prime},x^{\star})]$
for every true model $x^{\star}$. Accordingly, it is natural that the function
$\Psi$ also satisfies the property that
$\Psi(u,v)=\Psi(u^{\prime},v^{\prime})$; stated differently, one need only
specify $\Psi$ for a ‘minimal’ set of covering pairs. We present next a
definition that formalizes this notion precisely.

######  Definition 3 (Minimal covering pairs).

Consider a graded poset $(\mathcal{L},\preceq,\mathrm{rank}(\cdot))$ endowed
with a similarity valuation $\rho$. A subset
$\mathcal{S}\subset\\{(a,b)~{}|~{}b\text{ covers }a\text{ in }\mathcal{L}\\}$
of covering pairs in $\mathcal{L}$ is called _minimal_ if the following two
properties hold:

  * •

For each covering pair $(u^{\prime},v^{\prime})\notin\mathcal{S}$, there
exists $(u,v)\in\mathcal{S}$ with
$\mathrm{rank}(v)\leq\mathrm{rank}(v^{\prime})$ such that
$\rho(v,z)-\rho(u,z)=\rho(v^{\prime},z)-\rho(u^{\prime},z)$ for all
$z\in\mathcal{L}$.

  * •

For distinct covering pairs $(u,v),(u^{\prime},v^{\prime})\in\mathcal{S}$,
there exists some $z\in\mathcal{L}$ such that
$\rho(v,z)-\rho(u,z)\neq\rho(v^{\prime},z)-\rho(u^{\prime},z)$.

In words, a minimal set of covering pairs $\mathcal{S}$ for a graded poset
$\mathcal{L}$ is an inclusion-minimal collection of smallest rank covering
pairs for which it suffices to consider the values of $\Psi$. For Example 1 on
variable selection with the similarity valuation $\rho_{\text{meet}}$, a
minimal set of covering pairs is given by
$\mathcal{S}=\\{(\emptyset,\\{i\\})~{}|~{}i=1,\dots,p\\}$ and this minimal set
is unique. In general, however, such sets are not unique; see SI Appendix,
Section 3 where we derive minimal sets of covering pairs for several examples.
Minimal sets of covering pairs are significant methodologically from both
computational and statistical perspectives. In particular, several of our
bounds for discrete posets depend on the cardinality $|\mathcal{S}|$ and these
also involve computations that scale in number of operations with
$|\mathcal{S}|$. Therefore, identifying a minimal set of covering pairs that
is small in cardinality is central to the success of our proposed methods. In
the remainder of this section, we assume that a minimal set of covering pairs
$\mathcal{S}$ for a given model poset $\mathcal{L}$ is available.

###  3.2 Model Selection via Stability

Our first method for designing a suitable function $\Psi$ to employ in
Algorithm 1 is based on subsampling and corresponding model averaging. We
assume that we have access to a base procedure $\hat{x}_{\mathrm{base}}$ that
provides model estimates from data as well as a dataset $\mathcal{D}$
consisting of observations drawn from a probability distribution parameterized
by the true model $x^{\star}$, and our approach is to aggregate the model
estimates provided by $\hat{x}_{\mathrm{base}}$ on subsamples of
$\mathcal{D}$. The requirements on the quality of the procedure
$\hat{x}_{\mathrm{base}}$ are quite mild, and we prove bounds in the sequel on
the false discovery associated with the aggregated model. In particular, the
aggregation method ensures that the averaged model is ‘stable’ in the sense
that it contains discoveries that are supported by a large fraction of the
subsamples. Our method generalizes the stability selection method for variable
selection [14, 23] and subspace stability selection for subspace estimation
[25]. We demonstrate the broad applicability of this methodology in Section 4
by applying it to several examples from Section 2.

Formally, fix a positive even integer $B$ and obtain $B/2$ complementary
partitions of the dataset $\mathcal{D}$, each of which partitions
$\mathcal{D}$ into two subsamples of equal size. Let this collection of
subsamples be denoted $\\{\mathcal{D}^{(\ell)}\\}_{\ell=1}^{B}$, and let
$\hat{x}_{\mathrm{base}}(\mathcal{D}^{(\ell)})$ denote the model estimate
obtained by applying the base procedure to the subsample
$\mathcal{D}^{(\ell)}$. For any covering pair $(u,v)$ of a model poset
$\mathcal{L}$, we define:

$\Psi_{\mathrm{stable}}(u,v):=1-\frac{1}{B}\sum_{\ell=1}^{B}\frac{\rho(v,\hat{x}_{\mathrm{base}}(\mathcal{D}^{(\ell)}))-\rho(u,\hat{x}_{\mathrm{base}}(\mathcal{D}^{(\ell)}))}{c_{\mathcal{L}}(u,v)},$
(6)

where $c_{\mathcal{L}}(u,v):=\max_{z\in\mathcal{L}}\rho(v,z)-\rho(u,z)$.
Appealing to properties of similarity valuations, we have that
$\rho(v,\hat{x}_{\mathrm{base}}(\mathcal{D}^{(\ell)}))-\rho(u,\hat{x}_{\mathrm{base}}(\mathcal{D}^{(\ell)}))\geq
0$ and $c_{\mathcal{L}}(u,v)\geq 1$. The term
$\rho(v,\hat{x}_{\mathrm{base}}(\mathcal{D}^{(\ell)}))-\rho(u,\hat{x}_{\mathrm{base}}(\mathcal{D}^{(\ell)}))$
measures the additional discovery about
$\hat{x}_{\mathrm{base}}(\mathcal{D}^{(\ell)})$ in the model $v$ relative to
the model $u$, while the quantity $c_{\mathcal{L}}(u,v)$ serves as
normalization to ensure that $\Psi_{\mathrm{stable}}(u,v)\in[0,1]$. In
particular, $\Psi_{\mathrm{stable}}(u,v)$ being small implies that the
additional discovery represented by the model $v$ over the model $u$ is
supported by a large fraction of the subsamples
$\\{\mathcal{D}^{(\ell)}\\}_{\ell=1}^{B}$. Consequently, when
$\Psi_{\mathrm{stable}}$ is employed in the context of Algorithm 1 in which we
greedily grow a path, each ‘step’ in the path corresponds to a discovery that
is supported by a large fraction of the subsamples. We provide theoretical
support for this approach in Theorem 10 in the sequel and the proof proceeds
by showing that $\Psi_{\mathrm{stable}}(u,v)$ being small implies that
$\mathbb{E}[1-(\rho(u,x^{\star})-\rho(v,x^{\star}))]$ is small; we combine
this observation with the telescoping sum formula (5) to obtain a bound on the
expected false discovery of the model estimated by Algorithm 1.

When Algorithm 1 with $\Psi=\Psi_{\mathrm{stable}}$ is specialized to Example
1 and Example 8, we obtain the stability selection procedure of [14, 23] and
the subspace stability selection method of [25]. For variable selection in
particular, Algorithm 1 with $\Psi=\Psi_{\mathrm{stable}}$ outputs the subset
of variables that appear in at least a $1-\alpha$ fraction of the models
estimated by the base procedure when applied to the subsamples
$\\{\mathcal{D}^{(\ell)}\\}_{\ell=1}^{B}$. More generally, Algorithm 1 with
$\Psi=\Psi_{\mathrm{stable}}$ also provides a procedure for model selection in
Examples 2-7 corresponding to discrete model posets.

######  Theorem 10 (false discovery control for Algorithm 1 with
$\Psi=\Psi_{\mathrm{stable}}$).

Let $(\mathcal{L},\preceq,\mathrm{rank}(\cdot))$ be a graded discrete model
poset with integer-valued similarity valuation $\rho$ and let $\mathcal{S}$ be
an associated set of minimal covering pairs. Let $\hat{x}_{\mathrm{base}}$ be
a base estimator. Suppose the dataset $\mathcal{D}$ employed in the
computation of $\Psi_{\text{stable}}$ consists of i.i.d. observations from a
distribution parametrized by the true model $x^{\star}\in\mathcal{L}$, and
suppose $\hat{x}_{\mathrm{sub}}$ is an estimator obtained by applying
$\hat{x}_{\mathrm{base}}$ to a subsample of $\mathcal{D}$ of size
$|\mathcal{D}|/2$. Fix $\alpha\in(0,1/2)$ and a positive, even integer $B$.
The output $\hat{x}_{\text{stable}}$ from Algorithm 1 with
$\Psi=\Psi_{\mathrm{stable}}$ satisfies the following false discovery bound

$\mathbb{E}[\mathrm{FD}(\hat{x}_{\mathrm{stable}},x^{\star})]\leq\sum_{\begin{subarray}{c}(u,v)\in\mathcal{S}\cap\mathcal{T}_{\text{null}}\end{subarray}}\frac{\mathbb{E}[\rho(v,\hat{x}_{\mathrm{sub}})-\rho(u,\hat{x}_{\mathrm{sub}})]^{2}}{(1-2\alpha)c_{\mathcal{L}}(u,v)^{2}}.$
(7)

Here the set $\mathcal{T}_{\text{null}}:=\\{(u,v)\text{ covering pair in
}\mathcal{L}~{}|~{}\rho(v,x^{\star})=\rho(u,x^{\star})\\}$ consists of all
covering pairs $(u,v)$ for which there is no additional discovery in the model
$v$ over the model $u$ with respect to the true model $x^{\star}$.

In the bound (7), the numerator
$\mathbb{E}[\rho(v,\hat{x}_{\mathrm{sub}})-\rho(u,\hat{x}_{\mathrm{sub}})]^{2}$
of each summand characterizes the quality of the base estimator on subsamples;
base estimators for which this term is small, when employed in the computation
of $\Psi_{\text{stable}}$ in Algorithm 1, yield models
$\hat{x}_{\text{stable}}$ with small false discovery.

######  Remark 3.

When specialized to Example 1 on variable selection with similarity valuation
$\rho_{\text{meet}}$, we recover Theorem 1 of [23]. Specifically, in (7), we
have that $c_{\mathcal{L}}(u,v)=1$ for any covering pair $(u,v)$ and
$\sum_{(u,v)\in\mathcal{S}\cap\mathcal{T}_{\text{null}}}\mathbb{E}[\rho(v,\hat{x}_{\mathrm{sub}})-\rho(u,\hat{x}_{\mathrm{sub}})]^{2}=\sum_{\mathrm{null}~{}i}\mathbb{P}[\text{variable
}i\text{ selected by }\hat{x}_{\mathrm{sub}}]^{2}$.

Theorem 10 is general in its applicability to all the discrete posets in
Section 2, and it provides an intuitive bound on expected false discovery.
Nonetheless, it requires a characterization of the quality of the base
estimator $\hat{x}_{\mathrm{base}}$ employed on subsamples. When such a
characterization is unavailable, the false discovery bound (7) may not be
easily computable in practice. To address this shortcoming and obtain easily
computable bounds on false discovery, we consider natural assumptions on the
estimator $\hat{x}_{\mathrm{sub}}$ corresponding to the base estimator
$\hat{x}_{\mathrm{base}}$ applied to subsamples; these assumptions generalize
those developed in [14, 25] for stability-based methods for variable selection
and subspace estimation. To formulate these assumptions, we introduce some
notation. Let
$\mathrm{rank}(\mathcal{L}):=\max_{u\in\mathcal{L}}\mathrm{rank}(u)$ be the
largest rank of an element in $\mathcal{L}$ and let
$\mathcal{S}_{k}:=\\{(u,v)\in\mathcal{S}~{}|~{}\mathrm{rank}(v)=k\\}$ for each
$k\in[\mathrm{rank}(\mathcal{L})]$.

######  Assumption 1 (better than random guessing).

For each $k\in[\mathrm{rank}(\mathcal{L})]$ with
$\mathcal{S}_{k}\neq\emptyset$, we have that

$\displaystyle\sum_{\begin{subarray}{c}(u,v)\in\mathcal{S}_{k}\cap\mathcal{T}_{\mathrm{null}}\end{subarray}}\frac{1}{|\mathcal{S}_{k}\cap\mathcal{T}_{\mathrm{null}}|}\cdot\frac{\mathbb{E}[\rho(v,\hat{x}_{\mathrm{sub}})-\rho(u,\hat{x}_{\mathrm{sub}})]}{c_{\mathcal{L}}(u,v)}$
$\displaystyle\leq\sum_{\begin{subarray}{c}(u,v)\in\mathcal{S}_{k}\setminus\mathcal{T}_{\mathrm{null}}\end{subarray}}\frac{1}{|\mathcal{S}_{k}\setminus\mathcal{T}_{\mathrm{null}}|}\frac{\mathbb{E}[\rho(v,\hat{x}_{\mathrm{sub}})-\rho(u,\hat{x}_{\mathrm{sub}})]}{c_{\mathcal{L}}(u,v)}.$

######  Assumption 2 (invariance in mean).

For each $k\in[\mathrm{rank}(\mathcal{L})]$ with
$\mathcal{S}_{k}\neq\emptyset$, we have that
$\frac{\mathbb{E}[\rho(v,\hat{x}_{\mathrm{sub}})-\rho(u,\hat{x}_{\mathrm{sub}})]}{c_{\mathcal{L}}(u,v)}$
is the same for each $(u,v)\in\mathcal{S}_{k}\cap\mathcal{T}_{\mathrm{null}}$.

In words, Assumption 1 states that the average normalized difference in
similarity valuation of the estimator $\hat{x}_{\mathrm{sub}}$ is smaller over
‘null’ covering pairs than over non-null covering pairs. Assumption 2 states
that the expected value of the normalized difference in similarity of
$\hat{x}_{\mathrm{sub}}$ is the same for each ‘null’ covering pair. For the
case of variable selection (Example 1), Assumption 1 reduces precisely to the
‘better than random guessing’ assumption employed by [14], namely that the
expected number of true positives divided by the expected number of false
positives selected by the estimator $\hat{x}_{\mathrm{sub}}$ is larger than
the same ratio for an estimator that selects variables at random. As a second
condition, [14] required that the random variables in the collection
$\\{\mathbb{I}[i\in\hat{x}_{\mathrm{sub}}]:i\text{ null}\\}$ are exchangeable.
Our Assumption 2 when specialized to variable selection reduces to the
requirement that each of the random variables in the collection
$\\{\mathbb{I}[i\in\hat{x}_{\mathrm{sub}}]:i\text{ null}\\}$ has the same
mean. As a second illustration, consider the case of total ranking (Example 7)
involving items $a_{1},\dots,a_{p}$, with the least element
$\pi_{\text{null}}$ given by $\pi_{\text{null}}(a_{i})=i,~{}i=1,\dots,p$, the
true total ranking by $\pi^{\star}$, and the estimator on subsamples by
$\hat{\pi}_{\mathrm{sub}}$. Fix any $k\in\\{1,\dots,p-1\\}$. Assumption 1
states that the expected number of pairs
$(a_{i},a_{j})\in\text{inv}(\hat{\pi}_{\mathrm{sub}};\pi_{\text{null}})\cap\text{inv}(\pi^{\star};\pi_{\text{null}})$
with $j-i=k$ divided by the expected number of pairs
$(a_{i},a_{j})\in\text{inv}(\hat{\pi}_{\mathrm{sub}};\pi_{\text{null}})\setminus\text{inv}(\pi^{\star};\pi_{\text{null}})$
with $j-i=k$ is larger than the same ratio for an estimator that outputs a
total ranking at random. Assumption 2 states that the probability that
$(a_{i},a_{j})\in\text{inv}(\hat{\pi}_{\mathrm{sub}};\pi_{\text{null}})$ is
the same for all pairs $(a_{i},a_{j})$ with $j-i=k$ and
$(a_{i},a_{j})\not\in\text{inv}(\pi^{\star};\pi_{\text{null}})$. See SI
Appendix, Section 4 for a formal derivation.

######  Theorem 11 (refined false discovery control for Algorithm 1 with
$\Psi=\Psi_{\mathrm{stable}}$).

Consider the setup of Theorem 10, and suppose additionally that Assumptions 1
and 2 are satisfied. The output $\hat{x}_{\mathrm{stable}}$ from Algorithm 1
with $\Psi=\Psi_{\mathrm{stable}}$ satisfies the false discovery bound:

$\mathbb{E}[\mathrm{FD}(\hat{x}_{\mathrm{stable}},x^{\star})]\leq\sum_{k\in[\mathrm{rank}(\mathcal{L})],\mathcal{S}_{k}\neq\emptyset}\frac{q_{k}^{2}}{|\mathcal{S}_{k}|(1-2\alpha)},$
(8)

where
$q_{k}=\sum_{(u,v)\in\mathcal{S}_{k}}\mathbb{E}[\rho(v,\hat{x}_{\mathrm{sub}})-\rho(u,\hat{x}_{\mathrm{sub}})]/c_{\mathcal{L}}(u,v)$.

The quantities in the bound (8) may be readily computed in practice. In
particular, each $\mathcal{S}_{k}$ and $c_{\mathcal{L}}(\cdot,\cdot)$ depends
only on the model poset $\mathcal{L}$ and each $q_{k}$ can be approximated as
$q_{k}\approx\frac{1}{B}\sum_{\ell=1}^{B}\sum_{(u,v)\in\mathcal{S}_{k}}\frac{\rho(v,\hat{x}_{\mathrm{base}}(\mathcal{D}^{(\ell)}))-\rho(u,\hat{x}_{\mathrm{base}}(\mathcal{D}^{(\ell)}))}{c_{\mathcal{L}}(u,v)}$.
We give characterizations of the sets $\mathcal{S}_{k}$ and
$c_{\mathcal{L}}(\cdot,\cdot)$ for posets corresponding to total ranking,
partial ranking, clustering, and causal structure learning in SI Appendix,
Section 3.

######  Remark 4.

Specializing Theorem 11 to the case of variable selection, we arrive at the
bound in Theorem 1 of [14]. Specifically, note that for the Boolean poset with
the similarity valuation $\rho_{\text{meet}}$, $\mathcal{S}_{k}=\emptyset$ for
$k\geq 2$, $|\mathcal{S}_{1}|=\\#\text{ variables}$, and
$q_{1}=\sum_{i}\mathbb{P}[\text{variable }i\text{ selected by
}\hat{x}_{\mathrm{sub}}]$ is the average number of variables selected by the
estimator $\hat{x}_{\mathrm{sub}}$.

Turning our attention to Examples 8-9, the situation is considerably more
complicated with continuous model posets. A result for these two cases under
the same setup as in Theorem 10 yields the following bound for
$\alpha\in(0,1/2)$ (see SI Appendix, Section 5):

$\displaystyle\mathbb{E}[\mathrm{FD}(\hat{x}_{\mathrm{stable}},x^{\star})]~{}\leq{}$
$\displaystyle\frac{2\alpha+2\sqrt{\alpha}}{1-\alpha}\mathbb{E}[\mathrm{rank}(\hat{x}_{\mathrm{sub}})]+\mathbb{E}[\sqrt{\mathrm{FD}(\hat{x}_{\mathrm{sub}},x^{\star})}]^{2}.$
(9)

The first term in the bound is a function of the average number of discoveries
made by the estimator $\hat{x}_{\mathrm{sub}}$, and this term is smaller for
$\alpha\approx 0$. The second term in the bound concerns the quality of the
estimator $\hat{x}_{\mathrm{sub}}$. Specifically, note that Jensen’s
inequality implies
$\mathbb{E}[\sqrt{\mathrm{FD}(\hat{x}_{\mathrm{sub}},x^{\star})}]^{2}\leq\mathbb{E}[\mathrm{FD}(\hat{x}_{\mathrm{sub}},x^{\star})]$,
so that the improvement provided by the estimator $\hat{x}_{\mathrm{stable}}$
based on subsampling and model averaging over the estimator
$\hat{x}_{\mathrm{sub}}$ that simply employs the base estimator on subsamples
is characterized by
$\text{var}(\mathrm{FD}(\hat{x}_{\mathrm{sub}},x^{\star}))$. Thus, the key
remaining task as before is to characterize the properties of the estimator
$\hat{x}_{\mathrm{sub}}$. However, the difficulty with the continuous examples
is that conditions akin to Assumptions 1-2 are substantially more challenging
to formulate and analyze at an appropriate level of generality. (One such
effort under a limited setting for the case of subspace estimation is
described in [25].) It is of interest to develop such a general framework for
continuous model posets, and we leave this as a topic for future research.

###  3.3 Model Selection via Testing

Our second approach to designing a suitable function $\Psi$ to employ in
Algorithm 1 is based on testing the following null hypothesis for each
(minimal) covering pair $(u,v)$ of a discrete model poset $\mathcal{L}$:

$\displaystyle H_{0}^{u,v}$
$\displaystyle:\rho(v,x^{\star})=\rho(u,x^{\star}),$ (10)
$\displaystyle\Psi_{\mathrm{test}}(u,v)$ $\displaystyle:=\text{p-value
corresponding to }H_{0}^{u,v}.$

The null hypothesis $H_{0}^{u,v}$ in (10) states that there is no additional
discovery about $x^{\star}$ in the model $v$ relative to the model $u$, and
small values of $\Psi_{\mathrm{test}}(u,v)$ provide evidence for rejecting
this null hypothesis and accepting the alternative that
$\rho(v,x^{\star})>\rho(u,x^{\star})$. When $\Psi_{\mathrm{test}}$ is employed
in the context of Algorithm 1 in which we greedily grow a path, each ‘step’ in
the path corresponds to a discovery for which we have the ‘strongest evidence’
using the above test. Our next result provides theoretical support for this
method.

######  Theorem 12 (false discovery control for Algorithm 1 with
$\Psi=\Psi_{\mathrm{test}}$).

Let $(\mathcal{L},\preceq,\mathrm{rank}(\cdot))$ be a graded discrete model
poset with integer-valued similarity valuation $\rho$ and let $\mathcal{S}$ be
an associated set of minimal covering pairs. The output
$\hat{x}_{\mathrm{test}}$ of Algorithm 1 with $\Psi=\Psi_{\text{test}}$
satisfies the false discovery bound
$\mathbb{P}\left(\mathrm{FD}(\hat{x}_{\mathrm{test}},x^{\star})>0\right)\leq\alpha|\mathcal{S}|$.

The multiplicity factor involving the cardinality of the set of minimal
covering pairs $\mathcal{S}$ is akin to a Bonferroni-type correction, and it
highlights the significance of identifying a set of minimal covering pairs of
small cardinality. We emphasize that although Algorithm 1 with
$\Psi=\Psi_{\text{test}}$ proceeds via sequential hypothesis testing, the
procedure is applicable to general model classes with no underlying Boolean
logical structure; in particular, it is the graded poset structure underlying
our framework that facilitates such methodology.

As an illustration of the multiplicity factor $|\mathcal{S}|$ for different
settings, we have that $|\mathcal{S}|=p(p-1)$ for partial ranking;
$|\mathcal{S}|=\sum_{k=1}^{p-1}\big{(}\begin{smallmatrix}p\\\
k+1\end{smallmatrix}\big{)}\sum_{\ell=1}^{k}\big{(}\begin{smallmatrix}k+1\\\
\ell\end{smallmatrix}\big{)}$ for clustering; and
$|\mathcal{S}|=\frac{p(p-1)}{2}$ for total ranking. See SI Appendix, Section 3
for further details.

The graded poset structure of a model class can also yield more powerful model
selection procedures than those obtained by the greedy procedure of Algorithm
1. We give one such illustration next in which a collection of model estimates
that each exhibits zero false discovery (with high probability) can be
‘combined’ to derive a more complex model that also exhibits zero false
discovery. Formally, a poset $(\mathcal{L},\preceq)$ is said to possess a
_join_ if for each $x,y\in\mathcal{L}$ there exists a $z\in\mathcal{L}$
satisfying $(i)$ $z\succeq x,z\succeq y$ and $(ii)$ for any $w\in\mathcal{L}$
with $w\succeq x,w\succeq y$, we have $w\succeq z$; such a $z$ is called the
join of $x,y$ and posets that possess a join are called _join semi-lattices_
(these are dual to the notion of a meet defined in Section 2). Except for the
posets in Examples 4, 6, and 9, the posets in the other examples are join
semi-lattices (see SI Appendix, Section 1). For a model class that is a join
semi-lattice, suppose we are provided estimates
$\hat{x}^{(1)},\dots,\hat{x}^{(m)}$ of a true model $x^{\star}$ such that
$\mathrm{FD}(\hat{x}^{(j)},x^{\star})=0,~{}j=1,\dots,m$ (for example, by
appealing to greedy methods such as Algorithm 1 or its variants). Appealing to
the properties of a similarity valuation, we can conclude that the join
$\hat{x}_{\text{join}}$ of $\hat{x}^{(1)},\dots,\hat{x}^{(m)}$ satisfies
$\mathrm{FD}(\hat{x}_{\text{join}},x^{\star})=0$; in general,
$\mathrm{rank}(\hat{x}_{\text{join}})$ is larger than
$\mathrm{rank}(\hat{x}^{(1)}),\dots,\mathrm{rank}(\hat{x}^{(m)})$, and
therefore, this procedure is one way to obtain a more powerful model by
combining less powerful ones while still retaining control on the amount of
false discovery. The following result formalizes matters.

######  Proposition 13 (using joins to obtain more powerful models).

Let $(\mathcal{L},\preceq,\mathrm{rank}(\cdot))$ be a graded discrete model
poset that is a join semi-lattice with integer-valued similarity valuation
$\rho$ and let $\mathcal{S}$ be an associated set of minimal covering pairs.
Consider a collection of estimates $\hat{x}^{(1)},\dots,\hat{x}^{(m)}$ of a
true model $x^{\star}$ and let $\hat{x}_{\mathrm{join}}$ denote the join of
$\hat{x}^{(1)},\dots,\hat{x}^{(m)}$. Suppose for each
$\hat{x}^{(j)},~{}j=1,\dots,m$ there is a path from the least element of
$\mathcal{L}$ to $\hat{x}^{(j)}$ such that every covering pair $(u,v)$ along
the path satisfies $\Psi_{\mathrm{test}}(u,v)\leq\alpha$. Then we have the
false discovery bound
$\mathbb{P}(\mathrm{FD}(\hat{x}_{\mathrm{join}},x^{\star})>0)\leq\alpha|\mathcal{S}|$.

##  4 Experiments

We describe the results of numerical experiments on synthetic and real data in
this section. We employ Algorithm 1 with with both $\Psi=\Psi_{\text{stable}}$
and $\Psi=\Psi_{\text{test}}$. For the testing-based approach, the manner in
which p-values are obtained is described in the context of each application
and we set $\alpha$ equal to $0.05/|\mathcal{S}|$ for a given set
$\mathcal{S}$ of minimal covering pairs. For the stability-based approach, we
consider $B=100$ subsamples obtained by partitioning a given dataset $50$
times into subsamples of equal size and we set $\alpha=0.3$.

To obtain a desired level of expected false discovery with the stability-based
approach, we appeal to Theorem 11 as follows. In the bound (8), each $q_{k}$
can be derived by averaging over subsamples (as explained in the discussion
after the statement of Theorem 11) and all the other quantities are known. The
values of these $q_{k}$’s in turn depend on the model estimates returned by
the base procedure $\hat{x}_{\mathrm{base}}$ employed on the subsamples; in
particular, if the estimate is the least element then each $q_{k}$ equals
zero, and as $\hat{x}_{\mathrm{base}}$ returns models of increasing
complexity, the value of each $q_{k}$ generally increases. Building on this
observation, we tune parameters in $\hat{x}_{\mathrm{base}}$ to return
increasingly more complex models until the bound (8) is at the desired level.
For causal structure learning we employ Greedy Equivalence Search as our base
procedure with tuning via the regularization parameter that controls model
complexity [4]. For clustering, we employ $k$-means [11] as the base procedure
with tuning via the number of clusters. For our illustrations with ranking
problems (both partial and total) in which we are provided with pairwise
comparison data, our base procedure first employs the maximum-likelihood
estimator associated to the Bradley-Terry model [3], which returns a vector of
positive weights $\hat{w}$ of dimension equal to the number of items. Using
this $\hat{w}$ we associate numerical values to covering pairs; each covering
pair corresponds to increasing the complexity of a model by including a pair
of items $(i,j)$ to the inversion set (in total ranking) or to the relation
specifying a strict partial order (in partial ranking), and the value we
assign is the difference $\hat{w}_{j}-\hat{w}_{i}$. Our base procedure then
constructs a path starting from the least element by greedily adding covering
pairs of largest value at each step, provided these values are larger than a
regularization parameter $\lambda>0$; smaller values of $\lambda$ yield model
estimates of larger complexity, while larger values yield estimates of smaller
complexity.

Figure 2: Comparing the performance of Algorithm 1 with
$\Psi=\Psi_{\text{stable}}$ versus a non-subsampling approach for total
ranking, clustering, and causal structure learning. Each problem setting
corresponds to a pair of dots and a connecting line. The comparison is in
terms of the amount of false and true discoveries.

Finally, for causal structure learning, we restrict our search during the
model aggregation phase of Algorithm 1 to paths that yield CPDAG models in
which each connected component in the skeleton has diameter at most two; such
a restriction facilitates a simple characterization of covering pairs. This
restriction is not imposed on the output of the base procedure. Moreover, the
true model can be an arbitrary CPDAG.

###  4.1 Synthetic data

We describe experiments with synthetic data using Algorithm 1 with
$\Psi=\Psi_{\text{stable}}$.

Total ranking: We consider a total ranking problem with $p=30$ items. We
observe $n$ i.i.d. games between players $i,j$ with the outcome modeled as
$y_{ij\ell}\sim\texttt{Bernoulli}(w^{\star}_{i}/(w^{\star}_{i}+w^{\star}_{j}))$
for $\ell=1,\dots,n$, where $w^{\star}\in\mathbb{R}^{p}_{++}$ is a feature
vector and $n\in\\{200,250,300\\}$. We fix $w^{\star}$ by first defining
$\tilde{w}\in\mathbb{R}^{p}_{++}$ as $\tilde{w}_{i}=\tau^{i-1},~{}i=1,\dots,p$
for $\tau\in\\{0.97,0.98,0.99\\}$, and then setting $w^{\star}$ equal to a
permutation of $\tilde{w}$ in which we swap the entries $1,3$, the entries
$8,10$, the entries $15,17$, the entries $20,22$, and the entries $25,27$.
Smaller values of $\tau$ correspond to better-distinguished items, and hence
to easier problem instances. The base procedure is tuned such that the
expected false discovery in (8) is at most three.

Clustering: We consider a clustering problem with $p=20$ variables. The true
partition consists of $12$ clusters with five variables in one cluster,
another five variables in a second cluster, and the remaining variables in
singleton clusters. The $p$ variables are independent two-dimensional
Gaussians. Each variable in cluster $i$ has mean $(\mu_{i},0)$ and covariance
$\tfrac{1}{4}I$; each $\mu_{i}=i/d$ for $d\in\\{3,3.5,4\\}$. Smaller values of
$d$ correspond to better-separated clusters, and hence to easier problem
instances. We are provided $n$ i.i.d. observations of these variables for
$n\in\\{40,65,90\\}$. The base procedure is tuned such that the expected false
discovery in (8) is at most three.

Causal structure learning: We consider a causal structure learning problem
over $p=10$ variables. The true DAG is generated by considering a random total
ordering of the variables, drawing directed edges from higher nodes in the
ordering to lower nodes independently with probability $v\in\\{0.13,0.18\\}$,
and defining a linear structural causal model in which each variable is a
linear combination of its parents plus independent Gaussian noise with mean
zero and variance $\tfrac{1}{4}$. The coefficients in the linear combination
are drawn uniformly at random from the interval $[0.5,0.7]$. Larger values of
$v$ lead to denser DAGs, and hence to harder problem instances. We obtain $n$
i.i.d. observations from these models for
$n\in\\{1000,1200,1400,1600,1800\\}$. The base procedure is tuned such that
the expected false discovery in (8) is at most two.

For the preceding three problem classes, we compare the performance of our
stability-based methodology versus that of a non-subsampled approach in which
the base procedure (with suitable regularization) is applied to the entire
dataset. For total ranking, the non-subsampled procedure simply extracts the
ranking implied by the maximum-likelihood estimator associated to the Bradley-
Terry model. For clustering, the non-subsampled approach employs $k$-means
where the number of clusters is chosen to maximize the average silhoutte score
[20]. For causal structure learning, the non-subsampled approach applies
Greedy Equivalence Search with a regularization parameter chosen based on
holdout validation ($70\%$ of the data is used for training and the remaining
$30\%$ for validation). Figure 2 presents the results of our experiments
averaged over $50$ trials, and as the plots demonstrate, our stability-based
methods yield models with smaller false discovery than the corresponding non-
subsampled approaches. This reduction in false discovery comes at the expense
of a loss in power, which is especially significant for some of the harder
problem settings. However, in all cases our stability-based method provides
the desired level of control on expected false discovery.

###  4.2 Real data

We describe next experiments with real data.

Partial ranking of tennis players: We consider the task of partially ranking
six professional tennis players – Berdych, Djokovic, Federer, Murray, Nadal,
and Wawrinka – based on historical head-to-head matches of these players up to
the end of 2022. We apply Algorithm 1 with $\Psi=\Psi_{\text{stable}}$ and
with the base procedure tuned such that the expected false discovery in (8) is
at most three. The output of our procedure is a rank-nine model given by the
partial ranking {Djokovic, Nadal} $>$ {Berdych, Murray, Wawrinka} and
{Federer} $>$ {Berdych, Wawrinka}.

Total ranking of educational systems: We consider the task of totally ordering
$p=15$ OECD countries in reading comprehension based on test results from the
Programme for International Student Assessment (PISA). We take the null
ranking as the ordering of the countries based on performance in $2015$ (see
the first row in Table 2), and we wish to update this model based on 2018 test
scores (data obtained from [27]), with the number of test scores ranging from
$696$ to $3414$. We apply Algorithm 1 with $\Psi=\Psi_{\text{test}}$ and we
obtain p-values by modeling the average test score of each country as a
Gaussian. We set $\alpha=0.05/\tfrac{p(p-1)}{2}$ (here $\tfrac{p(p-1)}{2}$ is
the cardinality of a set of minimal covering pairs), which yields the
guarantee from Theorem 12 that the estimated model has zero false discovery
with probability at least $0.95$. The output of our procedure is the rank-nine
model given by the total ranking in the second row in Table 2.

| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
2015 base ranking | CAN | FIN | IRL | EST | KOR | JPN | NOR | NZL | DEU | POL | SvN | NLD | AUS | SWE | DNK
testing approach | FIN | IRL | EST | CAN | KOR | JPN | NOR | NZL | POL | DEU | AUS | SWE | SvN | DNK | NLD

Table 2: Ranking of nations according to PISA reading comprehension scores;
the first column is the 2015 ranking of $15$ OECD countries which serves as
the base ranking for our analysis: based on test results in 2018, we update
this ranking using Algorithm 1 based on $\Psi=\Psi_{\text{test}}$ with the
result shown in the second column.

Learning causal structure among proteins: We aim to learn causal relations
underlying $p=11$ phosphoproteins and phospholipids from a mass spectroscopy
dataset containing $854$ measurements of abundance levels in an observational
setting [21]. We apply Algorithm 1 with $\Psi=\Psi_{\mathrm{stable}}$ and with
the base procedure tuned such that the expected false discovery in (8) is at
most two. Figure 3 presents the rank-six CPDAG model obtained from our
algorithm and compares to the estimates obtained from the literature [15, 21,
28]. Our CPDAG estimate has fewer edges than those in [15, 21, 28], which do
not explicitly provide control on false discovery.

p38ERKMEKJNKPIP3PKAPKCPIP2PLcGAKTRAF

Edge |  [21]a |  [21]b | [16] | [6] |  [15]a |  [15]b | [28]
---|---|---|---|---|---|---|---
MEK – RAF | $\leftarrow$ | $\leftarrow$ | $\rightarrow$ | $\rightarrow$ |  | – | –
AKT – ERK |  | $\leftarrow$ | $\rightarrow$ | $\leftarrow$ | – | – | –
PIP3 – PIP2 | $\rightarrow$ | $\rightarrow$ | $\rightarrow$ | $\leftarrow$ | $\rightarrow$ | $\rightarrow$ | $\rightarrow$
AKT – PKA | $\leftarrow$ | $\leftarrow$ | $\leftarrow$ | $\leftarrow$ |  | $\leftarrow$ | $\rightarrow$
P38 $\rightarrow$ PKC | $\leftarrow$ | $\leftarrow$ | $\leftarrow$ | $\leftarrow$ |  | – |
JNK $\rightarrow$ PKC | $\leftarrow$ | $\leftarrow$ | $\leftarrow$ | $\leftarrow$ | $\leftarrow$ | – |

Figure 3: left: CPDAG obtained by Algorithm 1 with
$\Psi=\Psi_{\text{stable}}$; right: comparing the edges obtained by our
algorithm (shown in the leftmost column) with different causal discovery
methods (with indicated reference). The consensus network according to [21] is
denoted here by “[21]a” and their reconstructed network by “[21]b”; The
authors in [15] apply two methods, and the results are presented by “[15]a”
and “[15]b”. Here, “$-$” means that the edge direction is not identified.

##  5 Proofs

For notational ease, for a covering pair $(u,v)$ and element $z$ in the poset
$\mathcal{L}$, we define $f(u,v;z)\triangleq\rho(v,z)-\rho(u,z)$. Recall that
$\mathcal{T}_{\text{null}}\triangleq\\{(u,v)\text{ covering pair in
}\mathcal{L}~{}|~{}\rho(v,x^{\star})=\rho(u,x^{\star})\\}$. Our analysis
relies on the following lemmas with the proofs presented in Appendix F.

######  Lemma 14.

Fix a discrete model poset $\mathcal{L}$ with integer-valued similarity
valuation $\rho$. For any model $x\in\mathcal{L}$ with $(x_{0},\dots,x_{k})$
being any path from the least element $x_{0}=x_{\text{least}}$ to $x_{k}=x$,
we have that
$\mathrm{FD}(x,x^{\star})\leq\sum_{i=1}^{k}\mathbb{I}[(x_{i-1},x_{i})\in\mathcal{T}_{\text{null}}]$.
As a result, we have that $\mathrm{FD}(x,x^{\star})>0$ implies the existence
of some $i$ for which $(x_{i-1},x_{i})\in\mathcal{T}_{\text{null}}$.

######  Lemma 15.

For any covering pairs $(u,v)$ and $(x,y)$ with $v\preceq x$, we cannot have
that $f(u,v;z)=f(x,y;z)$ for all $z\in\mathcal{L}$.

###  5.1 Proof of Theorem 10

For notational convenience, we let
$\hat{x}^{(\ell)}_{\mathrm{base}}=\hat{x}_{\mathrm{base}}(\mathcal{D}^{(\ell)})$
where $\\{\mathcal{D}^{(\ell)}\\}_{\ell=1}^{B}$ are the subsamples of
$\mathcal{D}$. Let $\hat{x}_{\text{stable}}$ be the output of Algorithm 1 with
$\mathrm{rank}(\hat{x}_{\text{stable}})=\hat{k}$, and let
$(x_{0},\dots,x_{\hat{k}})$ be the associated path from the least element
$x_{0}=x_{\text{least}}$ to $x_{\hat{k}}=\hat{x}_{\text{stable}}$; we have
that
$\frac{1}{B}\sum_{\ell=1}^{B}{f(x_{i-1},x_{i};\hat{x}_{\mathrm{base}}^{(\ell)})}/{c_{\mathcal{L}}(x_{i-1},x_{i})}\geq(1-\alpha)$
for each $i=1,\dots,\hat{k}$. Let
$\mathcal{C}\triangleq\\{(x_{i-1},x_{i})~{}|~{}i=1,\dots,\hat{k}\\}$. From
Lemma 14, we also have that
$\mathrm{FD}(\hat{x}_{\text{stable}},x^{\star})\leq\left|\mathcal{C}\cap\mathcal{T}_{\text{null}}\right|$.
Combining these observations, we conclude that
$\mathrm{FD}(\hat{x}_{\text{stable}},x^{\star})\leq\sum_{(u,v)\in\mathcal{C}\cap\mathcal{T}_{\text{null}}}\mathbb{I}\left[\frac{1}{B}\sum_{\ell=1}^{B}\frac{f(u,v;\hat{x}_{\mathrm{base}}^{(\ell)})}{c_{\mathcal{L}}(u,v)}\geq
1-\alpha\right]$. Next, we observe that for each covering pair in
$\mathcal{C}$ there exists a covering pair in the minimal set $\mathcal{S}$
with the values of $f$ and $c_{\mathcal{L}}$ remaining the same; moreover,
distinct covering pairs in $\mathcal{C}$ map to distinct covering pairs in
$\mathcal{S}$ from Lemma 15. Thus, we conclude that
$\mathrm{FD}(\hat{x}_{\text{stable}},x^{\star})\leq\sum_{{(u,v)\in\mathcal{S}\cap\mathcal{T}_{\text{null}}}}\mathbb{I}\left[\frac{1}{B}\sum_{\ell=1}^{B}\frac{f(u,v;\hat{x}_{\mathrm{base}}^{(\ell)})}{c_{\mathcal{L}}(u,v)}\geq
1-\alpha\right]$. We then have the following sequence of steps:

$\displaystyle\mathrm{FD}(\hat{x}_{\text{stable}},x^{\star})$
$\displaystyle\leq\sum_{\begin{subarray}{c}(u,v)\in\\\
\mathcal{S}\cap\mathcal{T}_{\text{null}}\end{subarray}}\mathbb{I}\left[\frac{1}{B/2}\sum_{\ell=1}^{B/2}\sum_{i\in\\{0,1\\}}\frac{f(u,v;\hat{x}_{\mathrm{base}}^{(2\ell-i)})}{c_{\mathcal{L}}(u,v)}\geq
2-2\alpha\right]$ (11) $\displaystyle\leq\sum_{\begin{subarray}{c}(u,v)\in\\\
\mathcal{S}\cap\mathcal{T}_{\text{null}}\end{subarray}}\mathbb{I}\left[\frac{1}{B/2}\sum_{\ell=1}^{B/2}\prod_{i\in\\{0,1\\}}\frac{f(u,v;\hat{x}_{\mathrm{base}}^{(2\ell-i)})}{c_{\mathcal{L}}(u,v)}\geq
1-2\alpha\right].$

The second inequality follows from $ab\geq a+b-1$ for $a,b\in[0,1]$, where we
set $a=f(u,v;\hat{x}_{\mathrm{base}}^{(2\ell-1)})/{c_{\mathcal{L}}(u,v)}$ and
$b={f(u,v;\hat{x}_{\mathrm{base}}^{(2\ell)})}/{c_{\mathcal{L}}(u,v)}$, and
note that $f(u,v;z)/c_{\mathcal{L}}(u,v)\in[0,1]$ for any $z\in\mathcal{L}$.
Taking expectations on both sides of the preceding inequality, we finally seek
a bound on
$\mathbb{P}\left[\frac{1}{B/2}\sum_{\ell=1}^{B/2}\prod_{i\in\\{0,1\\}}\frac{f(u,v;\hat{x}_{\mathrm{base}}^{(2\ell-i)})}{c_{\mathcal{L}}(u,v)}\geq
1-2\alpha\right]$. We have that:

$\displaystyle\mathbb{P}\left[\frac{1}{B/2}\sum_{\ell=1}^{B/2}\prod_{i\in\\{0,1\\}}\frac{f(u,v;\hat{x}_{\mathrm{base}}^{(2\ell-i)})}{c_{\mathcal{L}}(u,v)}\geq
1-2\alpha\right]$
$\displaystyle\leq\frac{\mathbb{E}\left[\frac{1}{B/2}\sum_{\ell=1}^{B/2}\prod_{i\in\\{0,1\\}}\frac{f(u,v;\hat{x}_{\mathrm{base}}^{(2\ell-i)})}{c_{\mathcal{L}}(u,v)}\right]}{1-2\alpha}$
(12)
$\displaystyle=\frac{\mathbb{E}\left[f(u,v;\hat{x}_{\mathrm{sub}})\right]^{2}}{c_{\mathcal{L}}(u,v)^{2}(1-2\alpha)}.$

Here $\hat{x}_{\mathrm{sub}}$ represents the estimator corresponding to the
base procedure $\hat{x}_{\mathrm{base}}$ applied to a subsample of
$\mathcal{D}$ of size $|\mathcal{D}|/2$. The inequality follows from Markov’s
inequality, and the equality follows by noting that complementary bags are
independent and identically distributed. Combining (11) and (12), we obtain
the desired result.

###  5.2 Proof of Theorem 11

We have from Theorem 10 that:

$\mathbb{E}[\mathrm{FD}(\hat{x}_{\text{stable}},x^{\star})]\leq\sum_{k=1}^{\mathrm{rank}(\mathcal{L})}\sum_{(u,v)\in\mathcal{S}_{k}\cap\mathcal{T}_{\text{null}}}\frac{\mathbb{E}[f(u,v;\hat{x}_{\mathrm{sub}})]^{2}}{(1-2\alpha)c_{\mathcal{L}}(u,v)^{2}}.$

Our goal is to bound
$\mathbb{E}[f(u,v;\hat{x}_{\mathrm{sub}})]/c_{\mathcal{L}}(u,v)$ for
$(u,v)\in\mathcal{S}_{k}\cap\mathcal{T}_{\text{null}}$. Note that each $q_{k}$
may be decomposed as

$\displaystyle q_{k}$ $\displaystyle=\sum_{\begin{subarray}{c}(u,v)\in\\\
\mathcal{S}_{k}\cap\mathcal{T}_{\text{null}}\end{subarray}}\frac{\mathbb{E}[f(u,v;\hat{x}_{\mathrm{sub}})]}{c_{\mathcal{L}}(u,v)}+\sum_{\begin{subarray}{c}(u,v)\in\\\
\mathcal{S}_{k}\setminus\mathcal{T}_{\text{null}}\end{subarray}}\frac{\mathbb{E}[f(u,v;\hat{x}_{\mathrm{sub}})]}{c_{\mathcal{L}}(u,v)}.$

Appealing to Assumption 1, we have that

$\displaystyle q_{k}$
$\displaystyle\geq\left(1+\frac{|\mathcal{S}_{k}\setminus\mathcal{T}_{\text{null}}|}{|\mathcal{S}_{k}\cap\mathcal{T}_{\text{null}}|}\right)\sum_{(u,v)\in\mathcal{S}_{k}\cap\mathcal{T}_{\text{null}}}\frac{\mathbb{E}[f(u,v;\hat{x}_{\mathrm{sub}})]}{c_{\mathcal{L}}(u,v)}.$

Rearranging the terms, we obtain that

$\sum_{(u,v)\in\mathcal{S}_{k}\cap\mathcal{T}_{\text{null}}}\frac{\mathbb{E}[f(u,v;\hat{x}_{\mathrm{sub}})]}{c_{\mathcal{L}}(u,v)}\leq\frac{q_{k}}{|\mathcal{S}_{k}|}|\mathcal{S}_{k}\cap\mathcal{T}_{\text{null}}|.$

Appealing to Assumption 2, we have for each
$(u,v)\in\mathcal{S}_{k}\cap\mathcal{T}_{\text{null}}$ that
$\frac{\mathbb{E}[f(u,v;\hat{x}_{\mathrm{sub}})]}{c_{\mathcal{L}}(u,v)}\leq\frac{q_{k}}{|\mathcal{S}_{k}|}$.
Plugging this bound into the conclusion of Theorem 10 yields the desired
result.

###  5.3 Proof of Theorem 12

Let $\hat{x}_{\text{test}}$ be the output of Algorithm 1 with
$\mathrm{rank}(\hat{x}_{\text{test}})=\hat{k}$, and let
$(x_{0},\dots,x_{\hat{k}})$ be the associated path from the least element
$x_{0}=x_{\text{least}}$ to $x_{\hat{k}}=\hat{x}_{\text{test}}$; we have that
$\Psi_{\text{test}}(x_{i-1},x_{i})\leq\alpha$ for each $i=1,\dots,\hat{k}$.
Let $\mathcal{C}\triangleq\\{(x_{i-1},x_{i})~{}|~{}i=1,\dots,\hat{k}\\}$. From
Lemma 14, we have that $\mathrm{FD}(\hat{x}_{\text{test}},x^{\star})>0$
implies the existence of a covering pair
$(u,v)\in\mathcal{C}\cap\mathcal{T}_{\text{null}}$ for which
$\Psi_{\text{test}}(u,v)\leq\alpha$. For each covering pair in $\mathcal{C}$,
there exists a covering pair in $\mathcal{S}$ with the same value of
$\Psi_{\text{test}}$; thus, there exists
$(u,v)\in\mathcal{S}\cap\mathcal{T}_{\text{null}}$ such that
$\Psi_{\text{test}}(u,v)\leq\alpha$. Consequently:

$\displaystyle\mathbb{P}(\mathrm{FD}(\hat{x}_{\text{test}},x^{\star})>0)$
$\displaystyle\leq\mathbb{P}\left(\exists(u,v)\in\mathcal{S}\cap\mathcal{T}_{\text{null}}\text{
s.t. }\Psi_{\text{test}}(u,v)\leq\alpha\right)$ (13)
$\displaystyle\leq\sum_{(u,v)\in\mathcal{S}\cap\mathcal{T}_{\text{null}}}\mathbb{P}(\Psi_{\text{test}}(u,v)\leq\alpha)\leq\alpha|\mathcal{S}|.$

Here the second inequality follows from the union bound and the final
inequality follows from the fact that the random variable
$\Psi_{\text{test}}(u,v)$ is a valid p-value under the null hypothesis
$\rho(v,x^{\star})=\rho(u,x^{\star})$.

###  5.4 Proof of Proposition 13

For each $\hat{x}^{(j)},~{}j=1,\dots,m$, we are given that there is a path
from $x_{\text{least}}$ to $\hat{x}^{(j)}$ such that $\Psi_{\text{test}}$ is
bounded by $\alpha$ for each covering pair in the path; let
$\mathcal{C}^{(j)}$ be the set of these covering pairs. As described in
Section 3.3 in the discussion preceding Proposition 13,
$\mathrm{FD}(\hat{x}_{\text{join}},x^{\star})>0$ implies that
$\mathrm{FD}(\hat{x}^{(j)},x^{\star})>0$ for some $j=1,\dots,m$, which in turn
implies from Lemma 14 the existence of a covering pair
$(u,v)\in\mathcal{C}^{(j)}\cap\mathcal{T}_{\text{null}}$ for some
$j=1,\dots,m$. Following the same logic as in the proof of Theorem 12, we
conclude that $\mathrm{FD}(\hat{x}_{\text{join}},x^{\star})>0$ implies the
existence of $(u,v)\in\mathcal{S}\cap\mathcal{T}_{\text{null}}$ such that
$\Psi_{\text{test}}(u,v)\leq\alpha$. Using the same reasoning as in (13), we
have the desired conclusion.

##  6 Discussion

We present a general framework to endow a collection of models with poset
structure. This framework yields a systematic approach for quantifying model
complexity and false positive error in an array of complex model selection
tasks in which models are not characterized by Boolean logical structure (such
as in variable selection). Moreover, we develop methodology for controlling
false positive error in general model selection problems over posets, and we
describe experimental results that demonstrate the utility of our framework.

We finally discuss some future research questions that arise from our work. On
the mathematical front, a basic open question is to characterize fundamental
tradeoffs between false positive and false negative errors that are achievable
by any procedure in model selection over a general poset; this would
generalize the Neyman-Pearson lemma on optimal procedures for testing between
two hypotheses. On the computational and methodological front, it is of
interest to develop new methods to control false positive error as well as
false discovery rates, including in settings involving continuous model
posets.

## Acknowledgements

We thank Marina Meila and Lior Pachter for insightful conversations. AT
received funding from the Royalty Research Fund at the University of
Washington. PB received funding from the European Research Council under the
European Union’s Horizon 2020 research and innovation program (grant agreement
No. 786461). VC was supported in part by the Air Force Office of Scientific
Research grants FA9550-22-1-0225 and FA9550-23-1-0204 and by the National
Science Foundation grant DMS 2113724.

## References

  * [1] H. Akaike.  A new look at the statistical model identification.  IEEE Transactions on Automatic Control, 19:716–723, 1974.
  * [2] 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, 57:289–300, 1995.
  * [3] R. A. Bradley and M. E. Terry.  Rank analysis of incomplete block designs: I. the method of paired comparisons.  Biometrika, 39:324, 1952.
  * [4] D. M. Chickering.  Optimal structure identification with greedy search.  Journal of Machine Learning Research, 3:507–554, 2002.
  * [5] M. Drton and M. H. Maathuis.  Structure learning in graphical modeling.  Annual Review of Statistics and Its Application, 4(1):365–393, 2017\.
  * [6] D. Eaton and K. Murphy.  Exact bayesian structure learning from uncertain interventions.  In Proceedings of the 10th International Conference on Artificial Intelligence and Statistics (AISTATS), 2007.
  * [7] R. Fisher.  Statistical Methods for Research Workers.  Oliver $\&$ Boyd, 1932.
  * [8] D. P. Foster and R. A. Stine.  $\alpha$‐investing: a procedure for sequential control of expected false discoveries.  Journal of the Royal Statistical Society: Series B, 70, 2008.
  * [9] J. J. Goeman and U. R. Mansmann.  Multiple testing on the directed acyclic graph of gene ontology.  Bioinformatics, 24 4:537–44, 2008.
  * [10] H. Hotelling.  The generalization of student’s ratio.  Annals of Mathematical Statistics, 2:54–65, 1931.
  * [11] S. P. Lloyd.  Least squares quantization in pcm.  IEEE Trans. Inf. Theory, 28:129–136, 1982.
  * [12] G. Lorden.  Procedures for reacting to a change in distribution.  Annals of Mathematical Statistics, 42:1897–1908, 1971.
  * [13] N. Meinshausen.  Hierarchical testing of variable importance.  Biometrika, 95:265–278, 2008.
  * [14] N. Meinshausen and P. Bühlmann.  Stability selection.  Journal of the Royal Statistical Society: Series B, 72, 2010.
  * [15] N. Meinshausen, A. Hauser, J. Mooij, J. Peters, P. Versteeg, and P. Bühlmann.  Methods for causal inference from gene perturbation experiments and validation.  Proceeding of National Academy of Sciences, 113:7361–7368, 2016\.
  * [16] J. Mooij and T. Heskes.  Cyclic causal discovery from continuous equilibrium data.  In Proceedings of the 29th Annual Conference on Uncertainty in Artificial Intelligence, pages 431–439, 2013.
  * [17] J. Neyman and E. S. Pearson.  On the use and interpretation of certain test criteria for purposes of statistical inference part I.  Biometrika, pages 175–240, 1928.
  * [18] K. Popper.  Logik der forschung: Zur erkenntnistheorie der modernen naturwissenschaft.  Springer, 1935.
  * [19] P. R. Rosenbaum.  Testing hypotheses in order.  Biometrika, 95:248–252, 2008.
  * [20] P. J. Rousseeuw.  Silhouettes: a graphical aid to the interpretation and validation of cluster analysis.  Journal of Computational and Applied Mathematics, 20:53–65, 1987\.
  * [21] K. Sachs, O. Perez, D. Lauffenburger, and G. Nolan.  Causal protein-signaling networks derived from multiparameter single-cell data.  Science, 308:523–529, 2005.
  * [22] G. Schwarz.  Estimating the Dimension of a Model.  The Annals of Statistics, 6(2):461 – 464, 1978.
  * [23] R. Shah and R. Samworth.  Variable selection with error control: Another look at stability selection.  Journal of the Royal Statistical Society: Series B, 75, 2013.
  * [24] R. P. Stanley.  Enumerative Combinatorics, volume 1 of Cambridge Studies in Advanced Mathematics.  Cambridge University Press, 2nd edition, 2011.
  * [25] A. Taeb, P. Shah, and V. Chandrasekaran.  False discovery and its control in low-rank estimation.  Journal of the Royal Statistical Society: Series B, 92:997–1027, 2020.
  * [26] T. Verma and J. Pearl.  Equivalence and synthesis of causal models.  In Uncertainty in Artificial Intelligence, pages 255–270, 1990\.
  * [27] K. Wang.  Learning tower, 2021.  https://github.com/kevinwang09/learningtower.
  * [28] K. D. Yang, A. Katoff, and C. Uhler.  Characterizing and learning equivalence classes of causal dags under interventions.  In International Conference on Machine Learning, 2018.
  * [29] D. Yekutieli.  Hierarchical false discovery rate–controlling methodology.  Journal of the American Statistical Association, 103:309–316, 2008\.

## Appendix

##  Appendix A Meet Semi-lattice and Join Semi-lattice Properties and Posets
in Examples 1-9

The Boolean poset (Example 1), partition poset (Examples 2-3), integer poset
(Example 5), permutation poset (Example 7), and subspace poset (Example 8) are
all known in the literature to be lattices (and consequently meet-semi and
join semi-lattices); see [24].

We next show that for Examples 6 and 9 associated with partial ranking and
blind-source separation, the corresponding posets are also meet semi-lattices.
Consider the partial ranking setting in Example 6. Let $\mathcal{R}_{1}$ and
$\mathcal{R}_{2}$ be two relations that are irreflexive, asymmetric, and
transitive. Recalling that the partial ordering is based on inclusion, it is
clear that the relations
$\mathcal{R}=\\{(a,b):(a,b)\in\mathcal{R}_{1},(a,b)\in\mathcal{R}_{2}\\}$ is
the unique largest rank element in the partial ranking poset such that
$\mathcal{R}\preceq\mathcal{R}_{1}$ and $\mathcal{R}\preceq\mathcal{R}_{2}$.
Furthermore, for any $\tilde{\mathcal{R}}$ with
$\tilde{\mathcal{R}}\preceq\mathcal{R}_{1}$ and
$\tilde{\mathcal{R}}\preceq\mathcal{R}_{2}$, we clearly have that
$\tilde{\mathcal{R}}\preceq\mathcal{R}$. Consider the blind-source separation
setting in Example 9. Let $x_{1}$ and $x_{2}$ be two sets of linearly
independent subsets of unit norm vectors. Recalling that the partial ordering
in the associated poset is based on inclusion, it is clear that the set
$y=x_{1}\cap x_{2}$ is the unique largest rank element in the partial ranking
poset such that $y\preceq x_{1}$ and $y\preceq x_{2}$. Furthermore, for every
$z$ with $z\preceq x_{1}$ and $z\preceq x_{2}$, we have that $z\preceq y$.

We show that the poset corresponding to causal structure learning setting
(Example 4) is not meet semi-lattice or join semi-lattice. As a
counterexample, consider the CPDAGs $\mathcal{C}_{i}$ for $i=1,2,3,4$ shown in
Figure 4. Notice that $\mathcal{C}_{3}\preceq\mathcal{C}_{1}$,
$\mathcal{C}_{3}\preceq\mathcal{C}_{2}$,
$\mathcal{C}_{4}\preceq\mathcal{C}_{1}$, and
$\mathcal{C}_{4}\preceq\mathcal{C}_{2}$. Notice also that $\mathcal{C}_{3}$
and $\mathcal{C}_{4}$ are both CPDAGs with the largest rank that are smaller
(in a partial order sense) than $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$. We
thus can conclude that the poset is not meet semi-lattice. Similarly,
$\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are both CPDAGs with the smallest rank
that are larger (in a partial order sense) than $\mathcal{C}_{3}$ and
$\mathcal{C}_{4}$. We thus can conclude that the poset is not join semi-
lattice.

We next show that the poset for Example 6 is not join semi-lattice with a
simple counterexample. Consider as an example elements $x_{1}=\\{(1,2)\\}$ and
$x_{2}=\\{(2,1)\\}$. Note that there does not exist an element $z$ such that
$x_{1}\preceq z$ and $x_{2}\preceq z$. Thus, the poset is not join semi-
lattice.

Finally, we show that the poset corresponding to blind-source separation
(Example 9) is not join semi-lattice. Consider a collection of $p+1$ rank-1
elements in this poset, each element consisting of a single $p$ dimensional
vector. Then, evidently, there cannot exist an element $z$ consisting of a set
of vectors that contains all of the vectors in the rank-1 elements, while
satisfying the linear independence condition.

$1$$2$$3$ (a) $\mathcal{C}_{1}$

$1$$2$$3$ (b) $\mathcal{C}_{2}$

$1$$2$$3$ (c) $\mathcal{C}_{3}$

$1$$2$$3$ (d) $\mathcal{C}_{4}$

Figure 4: Four CPDAGs. Here, CPDAGs $\mathcal{C}_{3}$ and $\mathcal{C}_{4}$
are both largest complexity models that are smaller (in partial order sense)
than $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$. Similarly, CPDAGs
$\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are the smallest complexity models
that are larger (in a partial order sense) than $\mathcal{C}_{3}$ and
$\mathcal{C}_{4}$.

##  Appendix B Proof that $\mathrm{Eq.}$ (1) is a Similarity Valuation
Function

Recall that

$\rho_{\text{meet}}(x,y)=\max_{z\preceq x,z\preceq y}\mathrm{rank}(z).$ (14)

By definition, $\rho_{\text{meet}}(\cdot,\cdot)$ is a symmetric function. We
will now show that it satisfies the three properties in Definition 1 for any
pair of elements $x,y\in\mathcal{L}$. For the first property, we can conclude
$\rho_{\text{meet}}(x,y)\geq 0$ since by definition, the rank function returns
a non-negative integer for all the elements in the poset. Again, because of
the property of the rank function in a graded poset, a feasible $z$
(satisfying the constraints $z\preceq x$, $z\preceq y$) will necessarily have
$\mathrm{rank}(z)\leq\min\\{\mathrm{rank}(x),\mathrm{rank}(y)\\}$. For the
second property, consider any $w\in\mathcal{L}$ with $x\preceq w$. Note that:

$\rho_{\text{meet}}(w,y)=\max_{z\preceq w,z\preceq y}\mathrm{rank}(z).$ (15)

Then, any feasible $z$ in (14) is also feasible in (15) by the transitive
property of posets. Therefore,
$\rho_{\text{meet}}(x,y)\leq\rho_{\text{meet}}(w,y)$. For the third property,
first note that if $x\preceq y$, then $z=x$ is feasible in (14) and thus
$\rho_{\text{meet}}(x,y)\geq\mathrm{rank}(x)$. Since also
$\rho_{\text{meet}}(x,y)\leq\mathrm{rank}(x)$ by the second property of
similarity valuations, we have that
$\rho_{\text{meet}}(x,y)=\mathrm{rank}(x)$. Now suppose that
$\rho_{\text{meet}}(x,y)=\mathrm{rank}(x)$. By (14), we conclude that there
exists a feasible $z$ ($z\preceq x$, $z\preceq y$) such that
$\mathrm{rank}(z)=\mathrm{rank}(x)$. By the property of the rank function, we
have that if $\mathrm{rank}(z)=\mathrm{rank}(x)$ and $z\preceq x$, then $z=x$.
Since we have additionally that $z\preceq y$, we conclude that $x\preceq y$.

##  Appendix C Specializing Bound $\mathrm{Eq.}$ (8) for Different Problem
Settings

###  C.1 Partial Ranking

Let $S=\\{a_{1},a_{2},\dots,a_{p}\\}$ be the set of $p$ elements. We use the
similarity valuation $\rho:=\rho_{\text{meet}}$ in $\mathrm{Eq.}$ (1) of the
main paper.

####  C.1.1 Characterizing $\mathcal{S}$ for Partial Ranking

We construct a set $\mathcal{S}$ satisfying the properties in Definition 3 of
the main paper. Specifically, we let:

$\mathcal{S}=\\{(a_{i},a_{j}):i\neq j\\},$

with $|\mathcal{S}_{1}|=p(p-1)$ and $\mathcal{S}_{k}=\emptyset$ for every
$k\geq 2$.

We will show that set $\mathcal{S}$ as constructed above satisfies Definition
3. First, consider any covering pair
$({u^{\prime}},{v^{\prime}})\notin\mathcal{S}$. Here, ${u^{\prime}}$ and
${v^{\prime}}$ are relations and ${v^{\prime}}=u^{\prime}\cup(a_{i},a_{j})$
for some $i\neq j$. Then, for any $z\in\mathcal{L}$, it is easy to see that

$\rho({v^{\prime}},z)-\rho({u^{\prime}},z)=\mathbb{I}[(a_{i},a_{j})\in
z]=\rho(v,z)-\rho(u,z),$

where $v=\\{(a_{i},a_{j})\\}$ and $u=\emptyset$. Clearly,
$\mathrm{rank}(v)\leq\mathrm{rank}({v^{\prime}})$.

To show the second property, consider covering pairs
$(\\{(a_{i},a_{j})\\},\emptyset)\in\mathcal{S}$ and
$(\\{(a_{k},a_{l})\\},\emptyset)\in\mathcal{S}$. By construction of the set
$\mathcal{S}$, $(a_{i},a_{j})\neq(a_{k},a_{l})$. Let $z=\\{(a_{i},a_{j})\\}$.
Then, it is straightforward to see that
$\rho(\\{(a_{i},a_{j})\\},z)-\rho(\emptyset,z)=1$ but
$\rho(\\{(a_{k},a_{l})\\},z)-\rho(\emptyset,z)=0$.

####  C.1.2 Characterizing $c_{\mathcal{L}}(x,y)$ for Covering Pair $(x,y)$

Since for any $z$, $\rho(y,z)-\rho(x,z)=\mathbb{I}((a_{i},a_{j})\in z)$ for
some $(a_{i},a_{j})$. Thus, $c_{\mathcal{L}}(x,y)=1$.

####  C.1.3 Refined False Discovery Bound for Partial Ranking

Let $\hat{x}_{\text{stable}}$ be output of Algorithm 1 with
$\Psi=\Psi_{\text{stable}}$. Then:

$\mathbb{E}[\mathrm{FD}(\hat{x}_{\text{stable}},x^{\star})]\leq\frac{q_{1}^{2}}{(1-2\alpha)p(p-1)},$

where

$q_{1}=\sum_{i\neq j}\mathbb{I}[(a_{i},a_{j})\in\hat{x}_{\text{sub}}].$

Here, $\hat{x}_{\text{sub}}$ is the estimated partial ranking from supplying
$n/2$ samples to the base estimator. We can use the following data-driven
approximation for $q_{1}$:$q_{1}\approx\frac{1}{B}\sum_{\ell=1}^{B}\sum_{i\neq
j}\mathbb{I}[(a_{i},a_{j})\in\hat{x}_{\text{base}}(\mathcal{D}^{(\ell)})]$
with $\hat{x}_{\text{base}}(\mathcal{D}^{(\ell)}),l=1,2,\dots,B$ representing
the estimates from subsampling.

###  C.2 Total Ranking

Let $S=\\{a_{1},a_{2},\dots,a_{p}\\}$ be the set of $p$ elements. Let
$\pi_{\text{null}}(a_{i})=i$ for every $i=1,2,\dots,p$. We use the similarity
valuation $\rho:=\rho_{\text{total-ranking}}$ in $\mathrm{Eq.}$ (2) of the
main paper. As each element in the poset corresponds to a function $\pi:S\to
S$, we will use this functional notation throughout.

####  C.2.1 Characterizing $\mathcal{S}$ for Total Ranking

We construct a set $\mathcal{S}$ satisfying the properties in Definition 3 of
the main paper. Initialize $\mathcal{S}=\emptyset$. Then, for every relation
$(a_{i},a_{j})$ with $i<j$, we augment $\mathcal{S}$ as follows:

$\mathcal{S}=\mathcal{S}\cup(\pi_{1},\pi_{2}),$

where $\pi_{1},\pi_{2}$ are covering pairs. Here, $\pi_{2}$ is any rank $j-i$
element in the poset with the relation $(a_{i},a_{j})$ in its corresponding
inversion set. Furthermore, we let $\pi_{1}$ be a rank $j-i-1$ element that is
covered by $\pi_{2}$ and does not contain $(a_{i},a_{j})$ in its inversion
set. Recalling that
$\mathcal{S}_{k}=\\{(\pi_{1},\pi_{2})\in\mathcal{S},\mathrm{rank}(\pi_{2})=k\\}$,
we have that for every $k=1,2,\dots,p-1$

$|\mathcal{S}_{k}|=p-k.$

We will show that set $\mathcal{S}$ as constructed above satisfies Definition
3. First, consider any covering pair
$(\tilde{\pi}_{1},\tilde{\pi}_{2})\notin\mathcal{S}$. Then by definition, the
corresponding inversion sets are nested, i.e.
$\text{inv}(\tilde{\pi}_{2};\pi_{\text{null}})\supseteq\text{inv}(\tilde{\pi}_{1};\pi_{\text{null}})$
with the difference being a single relation. We will denote this relation by
$(a_{i},a_{j})$ with $j>i$. Consider the covering pair
$(\pi_{1},\pi_{2})\in\mathcal{S}$ where $(a_{i},a_{j})$ is in the inversion
set of $\pi_{2}$ but not in the inversion set of $\pi_{1}$. Then, for any
$\pi$, we have that

$\rho(\pi_{2},\pi)-\rho(\pi_{1},\pi)=\mathbb{I}((a_{i},a_{j})\in\text{inv}(\pi;\pi_{\text{null}}))=\rho(\tilde{\pi}_{2},\pi)-\rho(\tilde{\pi}_{1},\pi).$

Furthermore, it is straightforward to check that
$\mathrm{rank}(\tilde{\pi}_{2})\geq j-i=\mathrm{rank}(\pi_{2})$. We have thus
shown that $\mathcal{S}$ satisfies the first property in Definition 3.

To show the second property, consider covering pairs
$(\pi_{1},\pi_{2})\in\mathcal{S}$ where the difference between the two
inversion sets is the relation $(a_{i},a_{j})$. Let
$(\pi_{3},\pi_{4})\in\mathcal{S}$ where the difference between the two
inversion sets is the relation $(a_{k},a_{l})$. By construction of the set
$\mathcal{S}$, $(a_{i},a_{j})\neq(a_{k},a_{l})$. Let $\pi$ be a permutation
with $(a_{i},a_{j})$ in its inversion set. Then, as desired,

$\rho(\pi_{2},\pi)-\rho(\pi_{1},\pi)=\mathbb{I}((a_{i},a_{j})\in\text{inv}(\pi;\pi_{\text{null}}))\neq\rho({\pi}_{4},\pi)-\rho({\pi}_{3},\pi).$

####  C.2.2 Characterizing $c_{\mathcal{L}}(\pi_{1},\pi_{2})$ for Covering
Pair $(\pi_{1},\pi_{2})$

Since for any $\pi$,
$\rho(\pi_{2},\pi)-\rho(\pi_{1},\pi)=\mathbb{I}((a_{i},a_{j})\in\text{inv}(\pi;\pi_{\text{null}}))$
for some pair of elements $(a_{i},a_{j})$, then
$c_{\mathcal{L}}(\pi_{1},\pi_{2})=1$.

####  C.2.3 Refined False Discovery Bound for Total Ranking

Let $\hat{\pi}_{\text{stable}}$ be output of Algorithm 1 with
$\Psi=\Psi_{\text{stable}}$. Then:

$\mathbb{E}[\mathrm{FD}(\hat{\pi}_{\text{stable}},\pi^{\star})]\leq\sum_{k=1}^{p-1}\frac{q_{k}^{2}}{(1-2\alpha)(p-k)},$

where

$q_{k}=\sum_{(\pi_{1},\pi_{2})\in\mathcal{S}_{k}}\mathbb{E}[\rho(\pi_{2},\hat{\pi}_{\text{sub}})-\rho(\pi_{1},\hat{\pi}_{\text{sub}})]=\sum_{(i,j),j-i=k}\left[\mathbb{I}[(a_{i},a_{j})\in\text{inv}(\hat{\pi}_{\text{sub}};\pi_{\text{null}})\right].$

Here, $\hat{\pi}_{\text{sub}}$ represents ranking from supplying $n/2$ samples
to the base estimator. We can use the following data-driven approximation for
$q_{k}$:$q_{k}\approx\frac{1}{B}\sum_{(i,j),j-i=k}\sum_{\ell=1}^{B}\left[\mathbb{I}[(a_{i},a_{j})\in\text{inv}(\hat{\pi}_{\text{base}}(\mathcal{D}^{(\ell)});\pi_{\text{null}})]\right]$,
where $\hat{\pi}_{\text{base}}(\mathcal{D}^{(\ell)})$ represents the total
ranking obtained by supplying the base estimator on dataset
$\mathcal{D}^{(\ell)}$.

###  C.3 Clustering

We have a collection of $p$ items $\\{a_{1},a_{2},\dots,a_{p}\\}$ that we wish
to cluster. We let $x_{0}=\\{\\{a_{1}\\},\\{a_{2}\\},\dots,\\{a_{p}\\}\\}$ be
the least element. As described in the main paper, will use the similarity
valuation $\rho:=\rho_{\text{meet}}$ defined in $\mathrm{Eq.}$ (1) of the main
paper. Since the clustering poset is meet semi-lattice, $\rho$ computes the
rank of the meet of two elements; in this setting, the meet $x\wedge z$ of
$x=\\{G_{1},\dots,G_{q}\\}$ and $z=\\{\tilde{G}_{1},\dots,\tilde{G}_{s}\\}$ is

$x\wedge z=\\{G_{i}\cap\tilde{G}_{j}:G_{i}\cap\tilde{G}_{j}\neq\emptyset\\}.$

Subsequently, $\rho(x,z)=\mathrm{rank}(x\wedge y)$ is $p-\\#\text{ groups in
}x\wedge z$, which can be equivalently expressed as:

$\rho(x,z)=\sum_{i,j:|G_{i}\cap\tilde{G}_{j}|\neq\emptyset}|G_{i}\cap\tilde{G}_{j}|-1.$

For sets $G_{1},G_{2}\subseteq\\{1,2,\dots,p\\}$ with $G_{1}\cap
G_{2}=\emptyset$, we define:

$\mathcal{R}_{G_{1},G_{2}}:=\\{\\{a_{1}\\},\\{a_{2}\\},\dots,\\{a_{p}\\}\\}\setminus\\{\\{a_{i}\\}:a_{i}\in
G_{1}\cup G_{2}\\}.$

####  C.3.1 Characterizing $\mathcal{S}$ for Clustering

We construct a set $\mathcal{S}$ satisfying the properties in Definition 3.
Initialize $\mathcal{S}=\emptyset$. Then, for every $k=1,2,\dots,p-1$ and
pairs of groups of variables $G_{1}\subseteq\\{a_{1},\dots,a_{p}\\}$ and
$G_{2}\subseteq\\{a_{1},\dots,a_{p}\\}$ with $|G_{1}|+|G_{2}|=k+1$ and
$G_{1}\cap G_{2}=\emptyset$, we generate covering pairs $(x,y)$ with
$y=\\{G_{1}\cup G_{2},\mathcal{R}_{G_{1},G_{2}}\\}$ and
$x=\\{G_{1},G_{2},\mathcal{R}_{G_{1},G_{2}}\\}$, and let

$S=S\cup(x,y).$

Recalling that $\mathcal{S}_{k}=\\{(x,y)\in\mathcal{S},\mathrm{rank}(y)=k\\}$,
it is straightforward to check that for every $k=1,2,\dots,p-1$

$|\mathcal{S}_{k}|=\begin{pmatrix}p\\\
k+1\end{pmatrix}\sum_{\ell=1}^{k}\begin{pmatrix}k+1\\\ l\end{pmatrix}.$

Here, the terms $\big{(}\begin{smallmatrix}p\\\ k+1\end{smallmatrix}\big{)}$
counts the number of possible items in $G_{1}\cup G_{2}$ and the term
$\sum_{\ell=1}^{k+1}\big{(}\begin{smallmatrix}k+1\\\
l\end{smallmatrix}\big{)}$ counts the number of possible configurations of the
group $G_{2}$. We will show that the constructed set $\mathcal{S}$ satisfies
Definition 3 of the main paper. Our analysis is based on the following lemma.

######  Lemma 10.

Consider the covering pairs $(x,y)$ with $x=\\{G_{1},G_{2},\dots,G_{q}\\}$ and
$y=\\{G_{1}\cup G_{2},G_{3},\dots,G_{q}\\}$ where
$G_{i}\subseteq\\{1,2,\dots,px\\}$ and $G_{i}\cap G_{j}=\emptyset$ for every
$i\neq j$. Let $(\tilde{x},\tilde{y})$ be covering pairs with
$\tilde{y}=\\{G_{1}\cup G_{2},\mathcal{R}_{G_{1},G_{2}}\\}$ and
$\tilde{x}=\\{G_{1},G_{2},\mathcal{R}_{G_{1},G_{2}}\\}$. Then, for every
$z\in\mathcal{L}$, $\rho(y,z)-\rho(x,z)=\rho(\tilde{y},z)-\rho(\tilde{x},z)$.

###### Proof of Lemma 10.

Let $z=\\{\tilde{G}_{1},\dots,\tilde{G}_{s}\\}$ with
$\tilde{G}_{i}\subseteq\\{a_{1},a_{2},\dots,a_{p}\\}$ and
$\tilde{G}_{i}\cap\tilde{G}_{j}=\emptyset$ for every $i\neq j$. Then:

$\rho(y,z)=\sum_{j:(G_{1}\cup G_{2})\cap\tilde{G}_{j}\neq\emptyset}|(G_{1}\cup
G_{2})\cap\tilde{G}_{j}|-1+\sum_{i\geq
3,j:G_{i}\cap\tilde{G}_{j}\neq\emptyset}|G_{i}\cap\tilde{G}_{j}|-1,$

and

$\rho(x,z)=\sum_{j:G_{1}\cap\tilde{G}_{j}\neq\emptyset}|G_{1}\cap\tilde{G}_{j}|-1+\sum_{j:G_{2}\cap\tilde{G}_{j}\neq\emptyset}|G_{2}\cap\tilde{G}_{j}|-1+\sum_{i\geq
3,j:G_{i}\cap\tilde{G}_{j}\neq\emptyset}|G_{i}\cap\tilde{G}_{j}|-1.$

Since $\mathcal{R}_{G_{1},G_{2}}$ consists of groups of size one, we have
that:

$\rho(\tilde{y},z)=\sum_{j:(G_{1}\cup
G_{2})\cap\tilde{G}_{j}\neq\emptyset}|(G_{1}\cup G_{2})\cap\tilde{G}_{j}|-1,$

and

$\rho(\tilde{x},z)=\sum_{j:G_{1}\cap\tilde{G}_{j}\neq\emptyset}|G_{1}\cap\tilde{G}_{j}|-1+\sum_{j:G_{2}\cap\tilde{G}_{j}\neq\emptyset}|G_{2}\cap\tilde{G}_{j}|-1.$

We thus can see that
$\rho(y,z)-\rho(x,z)=\rho(\tilde{y},z)-\rho(\tilde{x},z)$. ∎

Showing $\mathcal{S}$ satisfies Definition 3 With Lemma 10 at hand, we show
that out constructed $\mathcal{S}$ satisfies Definition 3 of the main paper.
We start with the first property. Consider any
$(u^{\prime},v^{\prime})\subseteq\mathcal{L}$. Without loss of generality, we
take $v^{\prime}=\\{G_{1}\cup G_{2},G_{3},\dots,G_{q}\\}$ and
$u^{\prime}=\\{G_{1},G_{2},\dots,G_{q}\\}$. We let $v=\\{G_{1}\cup
G_{2},\mathcal{R}_{G_{1},G_{2}}\\}$ and
$u=\\{G_{1},G_{2},\mathcal{R}_{G_{1},G_{2}}\\}$. Then, according to Lemma 10,
we have that $\rho(v^{\prime},z)-\rho(u^{\prime},z)=\rho(v,z)-\rho(u,z)$.
Furthermore, since $\mathrm{rank}(x)=p-\\#\text{ groups in }x$, we have that
$\mathrm{rank}(v)\leq\mathrm{rank}(v^{\prime})$. Thus, the first property of
$\mathcal{S}$ is satisfied. We demonstrate the second property. Consider any
$(u,v)\in\mathcal{S}$ and $(u^{\prime},v^{\prime})\in\mathcal{S}$ that are
different. Let $u=\\{G_{1},G_{2},\mathcal{R}_{G_{1},G_{2}}\\}$ and
$v=\\{G_{1}\cup G_{2},\mathcal{R}_{G_{1},G_{2}}\\}$. Additionally, let
$u^{\prime}=\\{G_{1}^{\prime},G_{2}^{\prime},\mathcal{R}_{G_{1}^{\prime},G_{2}^{\prime}}\\}$
and $v^{\prime}=\\{G^{\prime}_{1}\cup
G^{\prime}_{2},\mathcal{R}_{G^{\prime}_{1},G^{\prime}_{2}}\\}$. Since the
covering pairs $(u,v)$ and $(u^{\prime},v^{\prime})$ are different, there must
exist two items $a_{i},a_{j}$ such that either $(a_{i},a_{j})$ are grouped
together in $v$ but are not together in $u$ or $(a_{i},a_{j})$ are grouped
together in $v^{\prime}$ but are not together in $u^{\prime}$. Let
$z=\\{\\{a_{i},a_{j}\\},\mathcal{R}_{\\{a_{i}\\},\\{a_{j}\\}}\\}$. Since
$\rho(v,z)-\rho(u,z)=\mathbb{I}[(a_{i},a_{j})\text{ grouped together in
}v\text{ but not in }u]$ and
$\rho(v^{\prime},z)-\rho(u^{\prime},z)=\mathbb{I}[(a_{i},a_{j})\text{ grouped
together in }v^{\prime}\text{ but not in }u^{\prime}]$, we have that
$\rho(v,z)-\rho(u,z)\neq\rho(v^{\prime},z)-\rho(u^{\prime},z)$.

####  C.3.2 Characterizing $c_{\mathcal{L}}(u,v)$ for Covering Pair $(u,v)$

######  Lemma 11.

Let $v=\\{G_{1}\cup{G}_{2},\mathcal{R}_{G_{1},G_{2}}\\}$ and
$u=\\{G_{1},{G}_{2},\mathcal{R}_{G_{1},G_{2}}\\}$ be a covering pair
$(u,v)\in\mathcal{S}$. Then, $c_{\mathcal{L}}(u,v)=\min\\{|G_{1}|,|G_{2}|\\}$.

###### Proof of Lemma 11.

Let $z=\\{\tilde{G}_{1},\dots,\tilde{G}_{q}\\}$. Then, from proof of Lemma 10,
we have that:

$\rho(v,z)-\rho(u,z)=\left[\sum_{j:(G_{1}\cup
G_{2})\cap\tilde{G}_{j}\neq\emptyset}|(G_{1}\cup
G_{2})\cap\tilde{G}_{j}|-1\right]-\left[\sum_{j:G_{1}\cap\tilde{G}_{j}\neq\emptyset}|G_{1}\cap\tilde{G}_{j}|-1\right]-\left[\sum_{j:G_{2}\cap\tilde{G}_{j}\neq\emptyset}|G_{2}\cap\tilde{G}_{j}|-1\right].$

Let $I_{1}:=\\{j:\tilde{G}_{j}\cap G_{1}\neq\emptyset\\}$ and
$I_{2}:=\\{j:\tilde{G}_{j}\cap G_{2}\neq\emptyset\\}$. Then,

$\rho(v,z)-\rho(u,z)=\left[\sum_{j\in I_{1}\cup I_{2}}|(G_{1}\cup
G_{2})\cap\tilde{G}_{j}|-1\right]-\left[\sum_{j\in{I}_{1}}|G_{1}\cap\tilde{G}_{j}|-1\right]-\left[\sum_{j\in
I_{2}}|G_{2}\cap\tilde{G}_{j}|-1\right].$

Simple manipulations yield:

$\rho(v,z)-\rho(u,z)=\left[\sum_{j\in I_{1}\cap I_{2}}|(G_{1}\cup
G_{2})\cap\tilde{G}_{j}|-1\right]-\left[\sum_{j\in{I}_{1}\cap
I_{2}}|G_{1}\cap\tilde{G}_{j}|-1\right]-\left[\sum_{j\in I_{1}\cap
I_{2}}|G_{2}\cap\tilde{G}_{j}|-1\right].$

Clearly, if $I_{1}\cap I_{2}=\emptyset$, then $\rho(v,z)-\rho(u,z)=0$. Suppose
$I_{1}\cap I_{2}\neq\emptyset$. Then,

$\rho(v,z)-\rho(u,z)=|I_{1}\cap I_{2}|+\left[\sum_{j\in I_{1}\cap
I_{2}}|(G_{1}\cup
G_{2})\cap\tilde{G}_{j}|-|G_{1}\cap\tilde{G}_{j}|-|G_{2}\cap\tilde{G}_{j}|\right]=|I_{1}\cap
I_{2}|.$

Notice that $|I_{1}\cap I_{2}|\leq\min\\{|G_{1}|,|G_{2}|\\}$. Then, the upper
bound can be achieved by for example setting
$z=\\{N,\\{\\{a_{1}\\},\\{a_{2}\\},\dots,\\{a_{p}\\}\setminus{N}\\}$ with
$N=\\{(a_{i},a_{j}):a_{i}\in G_{1},a_{j}\in G_{2}\\}$. ∎

####  C.3.3 Refined False Discovery Bound for Clustering

Let $\hat{x}_{\text{stable}}$ be output of Algorithm 1 with
$\Psi=\Psi_{\text{stable}}$. Then:

$\mathbb{E}[\mathrm{FD}(\hat{x}_{\text{stable}},x^{\star})]\leq\sum_{k=1}^{p-1}\frac{q_{k}^{2}}{(1-2\alpha)\begin{pmatrix}p\\\
k+1\end{pmatrix}\sum_{\ell=1}^{k}\begin{pmatrix}k+1\\\ l\end{pmatrix}},$

where,

$\displaystyle\begin{aligned}
q_{k}&=\sum_{(u,v)\in\mathcal{S}_{k}}\frac{\mathbb{E}[\rho(v,\hat{x}_{\text{sub}})-\rho(u,\hat{x}_{\text{sub}})]}{c(u,v)}\\\
&=\sum_{\begin{subarray}{c}G_{1}\subseteq\\{a_{1},\dots,a_{p}\\},G_{2}\subseteq\\{a_{1},\dots,a_{p}\\}\\\
G_{1}\cap
G_{2}=\emptyset;|G_{1}|+|G_{2}|=k+1\end{subarray}}\frac{\mathbb{E}[\\#\text{
groups }\hat{G}_{j}\text{ in }\hat{x}_{\text{sub}}\text{ satisfying
}\hat{G}_{j}\cap G_{1}\neq\emptyset\text{ and }\hat{G}_{j}\cap
G_{2}\neq\emptyset]}{\min\\{|G_{1}|,|G_{2}|\\}}.\end{aligned}$

Here, $\hat{x}_{\text{sub}}$ represents clustering from supplying $n/2$
samples to the base estimator. We will use the following data-driven
approximation to estimate $q_{k}$

$q_{k}\approx\frac{1}{B}\sum_{\begin{subarray}{c}G_{1}\subseteq\\{a_{1},\dots,a_{p}\\},G_{2}\subseteq\\{a_{1},\dots,a_{p}\\}\\\
G_{1}\cap
G_{2}=\emptyset;|G_{1}|+|G_{2}|=k+1\end{subarray}}\sum_{\ell=1}^{B}\frac{\\#\text{
groups }\hat{G}_{j}\text{ in
}\hat{x}_{\text{base}}(\mathcal{D}^{(\ell)})\text{ satisfying }\hat{G}_{j}\cap
G_{1}\neq\emptyset\text{ and }\hat{G}_{j}\cap
G_{2}\neq\emptyset]}{\min\\{|G_{1}|,|G_{2}|\\}},$

with $\hat{x}_{\text{base}}(\mathcal{D}^{(\ell)})$ represents the partition
obtained from supplying $\mathcal{D}^{(\ell)}$ to the base estimator.

###  C.4 Causal Structure Learning

Throughout, we consider covering pairs $(\mathcal{C}_{u},\mathcal{C}_{v})$
where each connected component in the skeletons of
$\mathcal{C}_{u},\mathcal{C}_{v}$ have a diameter at most two. We denote this
set by $\mathcal{T}$. Note that for any covering pair
$(\mathcal{C}_{u},\mathcal{C}_{v})\in\mathcal{T}$, $\mathcal{C}_{v}$ is a
polytree. Throughout, we will use the similarity valuation
$\rho:=\rho_{\text{meet}}$. Our analysis in this section will build on the
following result.

######  Lemma 12.

Let $\mathcal{C}_{u}$ and $\mathcal{C}_{v}$ be two CPDAGs that are polytrees
with $\mathcal{C}_{u}\preceq\mathcal{C}_{v}$. Then, the following statements
hold:

  1. (a)

for any pairs of nodes $\mathcal{E}$, the set of DAGs that result from
removing edges among pairs $\mathcal{E}$ in any DAG $\mathcal{G}_{v}$ form a
Markov equivalence class.

  2. (b)

for every DAG $\mathcal{G}_{v}\in\mathcal{C}_{v}$, there exists a DAG
$\mathcal{G}_{u}\in\mathcal{C}_{u}$ such that $\mathcal{G}_{u}$ is a directed
subgraph of $\mathcal{G}_{v}$.

###### Proof of Lemma 12.

We first prove part (a). By the polytree assumption, it follows that for any
DAG $\mathcal{G}_{v}$ in the CPDAG $\mathcal{C}_{v}$, removing the edges among
pairs in $\mathcal{E}$ does not create any v-structures, and removes the same
(potentially empty) v-structures. That means that the collection of DAGs
obtained by taking any DAG in $\mathcal{C}_{v}$ and removing the edges between
the pairs of nodes $\mathcal{E}$ will have the same skeleton and same
v-structures, and are thus in the same Markov equivalence class.

We next prove part (b). Let $(i,j)$ be the pair of nodes that are connected in
$\mathcal{C}_{v}$ but not in $\mathcal{C}_{u}$. Recall that
$\mathcal{C}_{u}\preceq\mathcal{C}_{v}$ implies there exists a DAG
$\mathcal{G}_{u}\in\mathcal{C}_{u}$ and a DAG
$\mathcal{G}_{v}\in\mathcal{C}_{v}$ where $\mathcal{G}_{u}$ is a subgraph of
$\mathcal{G}_{v}$, where $\mathcal{G}_{u}$ does not have the edge among pairs
$(i,j)$. Appealing to the result in part (a), we have that removing the edge
$(i,j)$ from any other DAG in $\mathcal{C}_{v}$ results in a DAG in the same
equivalence class, which is $\mathcal{C}_{u}$. ∎

####  C.4.1 Characterizing $\mathcal{S}$ for Causal Structure Learning

We construct the set $\mathcal{S}$ as follows. Initialize
$\mathcal{S}=\emptyset$. For every reference node, and $k=1,\dots,p-1$, let
$\mathcal{C}_{y}$ be a CPDAG generated with $k$ edges, where every edge is
between the reference node and another node; no other edges can be added
without violating the condition that the largest undirected path has size less
than or equal to two. A consequence of Lemma 12 is that there are $k$ CPDAGs
$\mathcal{C}_{x_{1}},\dots,\mathcal{C}_{x_{k}}$ that form a covering pair with
$\mathcal{C}_{y}$. We then let

$\mathcal{S}=\mathcal{S}\cup(\mathcal{C}_{x_{i}},\mathcal{C}_{y}),$

for every $i=1,2,\dots,k$. Recall that
$\mathcal{S}_{k}:=\\{(\mathcal{C}_{x},\mathcal{C}_{y})\in\mathcal{S},\mathrm{rank}(\mathcal{C}_{y})=k\\}$.
Then,

$|\mathcal{S}_{k}|=p\begin{pmatrix}p-1\\\
k\end{pmatrix}\sum_{i\in\\{0,2\dots,k\\}}\begin{pmatrix}k\\\ i\end{pmatrix}.$

The result above follows from noting that for every reference node and $k$
other nodes, there are
$\sum_{i\in\\{0,2\dots,k\\}}\big{(}\begin{smallmatrix}k\\\
i\end{smallmatrix}\big{)}$ possible CPDAGs that are polytrees can formed by
connecting the $k$ nodes to the reference node; the factor
$p\big{(}\begin{smallmatrix}p-1\\\ k\end{smallmatrix}\big{)}$ comes from $p$
total possible reference nodes and $\big{(}\begin{smallmatrix}p-1\\\
k\end{smallmatrix}\big{)}$ possible set of $k$ nodes to connect to the
reference node.

We will show that the constructed set $\mathcal{S}$ satisfies Definition 3 of
the main paper. Our analysis is based on the following lemma.

######  Lemma 13.

Let $\mathcal{C}_{\tilde{y}}$ be a CPDAG that contains $m$ disconnected
subgraphs (both directed and undirected). Let $\mathcal{C}_{\tilde{y}_{i}}$ be
each disconnected subgraph for $i=1,2,\dots,m$. Then, for any CPDAG
$\mathcal{C}_{z}$,

$\rho(\mathcal{C}_{\tilde{y}},\mathcal{C}_{z})=\sum_{i=1}^{m}\rho(\mathcal{C}_{\tilde{y}_{i}},\mathcal{C}_{z}).$

###### Proof.

We will first show that
$\rho(\mathcal{C}_{\tilde{y}},\mathcal{C}_{z})\leq\sum_{i=1}^{m}\rho(\mathcal{C}_{\tilde{y}_{i}},\mathcal{C}_{z})$.
Let
$\mathcal{C}_{\tilde{x}}\in\operatorname*{argmax}_{\mathcal{C}_{x}\preceq\mathcal{C}_{\tilde{y}},\mathcal{C}_{x}\preceq\mathcal{C}_{z}}\mathrm{rank}(\mathcal{C}_{x})$.
By definition, $\mathcal{C}_{x}\preceq\mathcal{C}_{\tilde{y}}$ if there is a
DAG $\mathcal{G}_{x}$ in $\mathcal{C}_{x}$ and a DAG $\mathcal{G}_{\tilde{y}}$
in $\mathcal{C}_{\tilde{y}}$ such that $\mathcal{G}_{x}$ is a subgraph of
$\mathcal{G}_{\tilde{y}}$. Since $\mathcal{G}_{\tilde{y}}$ has disconnected
components, so must $\mathcal{G}_{x}$. We let $\mathcal{C}_{\tilde{x}_{i}}$ be
the subgraphs of $\mathcal{C}_{\tilde{x}}$ where every subgraph
$\mathcal{C}_{\tilde{x}_{i}}$ only contains edges among nodes that are
connected (to other nodes) in the graph $\mathcal{C}_{\tilde{y}_{i}}$. By
construction, $\mathcal{C}_{\tilde{x}_{i}}\preceq\mathcal{C}_{\tilde{y}_{i}}$,
$\mathrm{rank}(\mathcal{C}_{\tilde{x}})=\sum_{i=1}^{m}\mathrm{rank}(\mathcal{C}_{\tilde{x}_{i}})$,
and $\mathcal{C}_{\tilde{x}_{i}}\preceq\mathcal{C}_{z}$. Thus,
$\mathrm{rank}(\mathcal{C}_{\tilde{x}_{i}})\leq\rho(\mathcal{C}_{\tilde{y}_{i}},\mathcal{C}_{z})$.
Then, we can conclude that

$\sum_{i=1}^{m}\rho(\mathcal{C}_{\tilde{y}_{i}},\mathcal{C}_{z})\geq\sum_{i=1}^{m}\mathrm{rank}(\mathcal{C}_{\tilde{x}_{i}})=\mathrm{rank}(\mathcal{C}_{\tilde{x}})=\rho(\mathcal{C}_{\tilde{y}},\mathcal{C}_{z}).$

Now we will show that
$\rho(\mathcal{C}_{\tilde{y}},\mathcal{C}_{z})\geq\sum_{i=1}^{m}\rho(\mathcal{C}_{\tilde{y}_{i}},\mathcal{C}_{z})$.
Let
$\mathcal{C}_{\tilde{x}_{i}}\in\operatorname*{argmax}_{\mathcal{C}_{x}\preceq\mathcal{C}_{\tilde{y}_{i}},\mathcal{C}_{x}\preceq\mathcal{C}_{z}}\mathrm{rank}(\mathcal{C}_{x})$.
Now form a CPDAG $\mathcal{C}_{\bar{y}}$ by combining all the disjoint graphs
$\mathcal{C}_{\tilde{x}_{i}}$ for every $i=1,2,\dots,m$ into one graph. Since
these graphs are disjoint (i.e. nodes that are connected in each graph are
distinct), we have that $\mathcal{C}_{\bar{y}}\preceq\mathcal{C}_{\tilde{y}}$
and $\mathcal{C}_{\bar{y}}\preceq\mathcal{C}_{z}$ and that
$\mathrm{rank}(\mathcal{C}_{\bar{y}})=\sum_{i=1}^{m}\mathrm{rank}(\mathcal{C}_{\tilde{x}_{i}})$.
So we conclude that

$\rho(\mathcal{C}_{\tilde{y}},\mathcal{C}_{z})\geq\mathrm{rank}(\mathcal{C}_{\bar{y}})=\sum_{i=1}^{m}\mathrm{rank}(\mathcal{C}_{\tilde{x}_{i}})=\sum_{i=1}^{m}\rho(\mathcal{C}_{\tilde{y}_{i}},\mathcal{C}_{z}).$

∎

Showing $\mathcal{S}$ satisfies Definition 3 For the first property, consider
covering pairs $(\mathcal{C}_{u^{\prime}},\mathcal{C}_{v^{\prime}})\in T$. Let
$(i,j)$ be the pair of nodes that are connected in $\mathcal{C}_{v^{\prime}}$
and are not connected in $\mathcal{C}_{u^{\prime}}$. Since every undirected
path in $\mathcal{C}_{v^{\prime}}$ has size at most $2$, then
$\mathcal{C}_{v^{\prime}}$ decouples into two disconnected CPDAGs
$\mathcal{C}_{v}$ and $\mathcal{C}_{1}$, where $\mathcal{C}_{v}$ only involves
nodes adjacent to $(i,j)$. Similarly, $\mathcal{C}_{u^{\prime}}$ decouples
into two disconnected CPDAGs $\mathcal{C}_{u}$ and $\mathcal{C}_{2}$, where
$\mathcal{C}_{2}=\mathcal{C}_{1}$ and $\mathcal{C}_{u}$ is covered by
$\mathcal{C}_{v}$. From Lemma 13, we have that for any CPDAG $\mathcal{C}_{z}$

$\rho(\mathcal{C}_{v^{\prime}},\mathcal{C}_{z})-\rho(\mathcal{C}_{u^{\prime}},\mathcal{C}_{z})=\rho(\mathcal{C}_{v},\mathcal{C}_{z})-\rho(\mathcal{C}_{u},\mathcal{C}_{z}).$

Notice that $(\mathcal{C}_{u},\mathcal{C}_{v})\in\mathcal{S}$. Furthermore,
since the number of edges (directed and undirected) in
$\mathcal{C}_{v^{\prime}}$ is larger than $\mathcal{C}_{v}$, we have that
$\mathrm{rank}(\mathcal{C}_{v})\leq\mathrm{rank}(\mathcal{C}_{v^{\prime}})$.

We next show the second property in Definition 3. Let
$(\mathcal{C}_{u},\mathcal{C}_{v})\in\mathcal{S}$ and
$(\mathcal{C}_{u^{\prime}},\mathcal{C}_{v^{\prime}})\in\mathcal{S}$. Our
objective is to show that
$\rho(\mathcal{C}_{v},\mathcal{C}_{z})-\rho(\mathcal{C}_{u},\mathcal{C}_{z})=\rho(\mathcal{C}_{v^{\prime}},\mathcal{C}_{z})-\rho(\mathcal{C}_{u^{\prime}},\mathcal{C}_{z})$
for all
$\mathcal{C}_{z}\Leftrightarrow\mathcal{C}_{u}=\mathcal{C}_{u^{\prime}}$ and
$\mathcal{C}_{v}=\mathcal{C}_{v^{\prime}}$. The direction $\leftarrow$
trivially holds, and hence we focus on the direction $\rightarrow$. We
consider multiple scenarios; throughout the extra edge that is present in
$\mathcal{C}_{v}$ and not in $\mathcal{C}_{u}$ is between the pair of nodes
$(i,j)$, and the extra edge that is present in $\mathcal{C}_{v^{\prime}}$ and
not in $\mathcal{C}_{u^{\prime}}$ is between the pair of nodes $(k,l)$.

  * (1)

Suppose that the nodes $(k,l)$ are not connected in $\mathcal{C}_{v}$. Letting
$\mathcal{C}_{z}$ be a CPDAG with only an edge between nodes $(k,l)$, we find
that
$\rho(\mathcal{C}_{v},\mathcal{C}_{z})-\rho(\mathcal{C}_{u},\mathcal{C}_{z})=0$
and
$\rho(\mathcal{C}_{v^{\prime}},\mathcal{C}_{z})-\rho(\mathcal{C}_{u^{\prime}},\mathcal{C}_{z})=1$.
So this scenario cannot occur.

  * (2)

Suppose there is an edge between pairs $(s,t)$ in $\mathcal{C}_{u^{\prime}}$
that is missing in $\mathcal{C}_{v}$ (and as a result in $\mathcal{C}_{u}$).
Construct CPDAG $\mathcal{C}_{z}$ with two edges, one between the pair $(i,j)$
and another between the pair $(s,t)$ with the property that
$\mathcal{C}_{z}\not\preceq\mathcal{C}_{v^{\prime}}$; this construction is
possible since
$(\mathcal{C}_{u^{\prime}},\mathcal{C}_{v^{\prime}})\in\mathcal{S}$, meaning
that if there is an edge between pair of nodes $(i,j)$ in
$\mathcal{C}_{v^{\prime}}$, this edge is incident to the edge between the pair
of nodes $(s,t)$. Then, it is evident that
$\rho(\mathcal{C}_{v},\mathcal{C}_{z})-\rho(\mathcal{C}_{u},\mathcal{C}_{z})=1$
but
$\rho(\mathcal{C}_{v^{\prime}},\mathcal{C}_{z})-\rho(\mathcal{C}_{u^{\prime}},\mathcal{C}_{z})=0$.
So this scenario cannot occur.

  * (3)

Suppose there is an edge between pairs $(s,t)$ in $\mathcal{C}_{u^{\prime}}$
that is missing in $\mathcal{C}_{u}$ but is not missing in $\mathcal{C}_{v}$.
Let $\mathcal{C}_{z}$ be a CPDAG only containing an edge between $(s,t)$. Then
it follows that
$\rho(\mathcal{C}_{v},\mathcal{C}_{z})-\rho(\mathcal{C}_{u},\mathcal{C}_{z})=1$
but
$\rho(\mathcal{C}_{v^{\prime}},\mathcal{C}_{z})-\rho(\mathcal{C}_{u^{\prime}},\mathcal{C}_{z})=0$.
So this scenario cannot occur.

From the impossibilities of scenarios 1-2, and noting that a similar argument
can be made by swapping $\mathcal{C}_{u^{\prime}}$ with $\mathcal{C}_{u}$, and
$\mathcal{C}_{v^{\prime}}$ with $\mathcal{C}_{v}$, we conclude that
$\mathcal{C}_{v},\mathcal{C}_{v^{\prime}}$ have edges between the same pairs
of nodes. Combining this result with the impossibility of scenario 3, we
conclude that $\mathcal{C}_{u},\mathcal{C}_{u^{\prime}}$ have edges between
the same pairs of nodes. We then continue with the final scenario.

  * (4)

Suppose that $\mathcal{C}_{v}$ and $\mathcal{C}_{v^{\prime}}$ are not
identical CPDAGs. Since both $\mathcal{C}_{v}$ and $\mathcal{C}_{v^{\prime}}$
have maximum undirected path length less than or equal to two, they both must
have the same reference node $i$ (where the other nodes are connected to).
Furthermore, since $\mathcal{C}_{v}$ and $\mathcal{C}_{v^{\prime}}$ have the
same skeleton and are different, they must have strictly more than one edge,
and they must have different v-structures. As a first sub-case, suppose
$\mathcal{C}_{v^{\prime}}$ have a v-structure $s\rightarrow i\leftarrow t$
that is not present in $\mathcal{C}_{v}$, so that $s\leftarrow i$ or $s-i$ in
$\mathcal{C}_{v}$. Then, let $\mathcal{C}_{z}$ be a CPDAG containing two edges
between the pairs $(i,j)$ and $(i,s)$ with
$\mathcal{C}_{z}\preceq\mathcal{C}_{v}$. By construction,
$\rho(\mathcal{C}_{v},\mathcal{C}_{z})-\rho(\mathcal{C}_{u},\mathcal{C}_{z})=1$
but
$\rho(\mathcal{C}_{v^{\prime}},\mathcal{C}_{z})-\rho(\mathcal{C}_{u^{\prime}},\mathcal{C}_{z})=0$.
Swapping $\mathcal{C}_{u^{\prime}}$ with $\mathcal{C}_{u}$, and
$\mathcal{C}_{v^{\prime}}$ with $\mathcal{C}_{v}$, and following similar
arguments, we arrive again at a contradiction if $\mathcal{C}_{v}$ has a
v-structure that is not present in $\mathcal{C}_{v^{\prime}}$.

From the impossibility of scenario 4, we conclude that $\mathcal{C}_{v}$ and
$\mathcal{C}_{v^{\prime}}$ have the same skeleton and v-structure and
consequently $\mathcal{C}_{v}=\mathcal{C}_{v^{\prime}}$. We thus have that
$\mathcal{C}_{u}\preceq\mathcal{C}_{v}$ and
$\mathcal{C}_{u^{\prime}}\preceq\mathcal{C}_{v}$. Furthermore, since
$\mathcal{C}_{u^{\prime}}$ and $\mathcal{C}_{u}$ have the same skeleton, both
are missing an edge between pair of nodes $(i,j)$ that is connected in
$\mathcal{C}_{v}$. Appealing to part a of Lemma 12, we conclude that
$\mathcal{C}_{u}=\mathcal{C}_{u^{\prime}}$.

####  C.4.2 Characterizing $c_{\mathcal{L}}(\mathcal{C}_{u},\mathcal{C}_{v})$
for Covering Pairs $(\mathcal{C}_{u},\mathcal{C}_{v})$

We have the following lemma.

######  Lemma 14.

Let $(\mathcal{C}_{u},\mathcal{C}_{v})$ be CPDAGs that are polytrees and form
a covering pair. Then, $c_{\mathcal{L}}(\mathcal{C}_{u},\mathcal{C}_{v})=1$.

###### Proof.

Let the pair of nodes $(i,j)$ be connected in $\mathcal{C}_{v}$ and not
connected in $\mathcal{C}_{u}$. Consider any CPDAG $\mathcal{C}_{z}$. Let
$\mathcal{C}_{\tilde{y}}\in\operatorname*{argmax}_{\mathcal{C}_{y}\preceq\mathcal{C}_{v},\mathcal{C}_{y}\preceq\mathcal{C}_{z}}\mathrm{rank}(\mathcal{C}_{y})$.
Since the CPDAG $\mathcal{C}_{v}$ is a polytree, so is the CPDAG
$\mathcal{C}_{\tilde{y}}$. Let $\mathcal{G}_{v}$ be any DAG in
$\mathcal{C}_{v}$. Then, by Lemma 12, there exists DAGs
$\mathcal{G}_{\tilde{y}}^{(1)}\in\mathcal{C}_{\tilde{y}}$ and
$\mathcal{G}_{u}\in\mathcal{C}_{u}$ such that $\mathcal{G}_{\tilde{y}}^{(1)}$
and $\mathcal{G}_{u}$ are both subgraphs of $\mathcal{G}_{v}$. Suppose we
remove an edge that may be present between the pair of nodes $(i,j)$ in
$\mathcal{G}_{\tilde{y}}^{(1)}$ and denote the resulting subgraph by
$\mathcal{G}_{x}^{(1)}$. By construction, $\mathcal{G}_{x}^{(1)}$ is also a
subgraph of $\mathcal{G}_{u}$. Since
$\mathcal{C}_{\tilde{y}}\preceq\mathcal{C}_{z}$, there exists a DAG
$\mathcal{G}_{\tilde{y}}^{(2)}\in\mathcal{C}_{\tilde{y}}$ and a DAG
$\mathcal{G}_{z}\in\mathcal{C}_{z}$ such that $\mathcal{G}_{\tilde{y}}^{(2)}$
is a subgraph of $\mathcal{G}_{z}$. Suppose again we remove an edge that may
be present between the pair of nodes $(i,j)$ in
$\mathcal{G}_{\tilde{y}}^{(2)}$ and denote the resulting subgraph by