# Biophysical models of cis-regulation as
interpretable neural networks

Ammar Tareen
Simons Center for Quantitative Biology
Cold Spring Harbor Laboratory
Cold Spring Harbor, NY 11724
<EMAIL_ADDRESS>
Justin B. Kinney
Simons Center for Quantitative Biology
Cold Spring Harbor Laboratory
Cold Spring Harbor, NY 11724
<EMAIL_ADDRESS>
###### Abstract

The adoption of deep learning techniques in genomics has been hindered by the
difficulty of mechanistically interpreting the models that these techniques
produce. In recent years, a variety of post-hoc attribution methods have been
proposed for addressing this neural network interpretability problem in the
context of gene regulation. Here we describe a complementary way of
approaching this problem. Our strategy is based on the observation that two
large classes of biophysical models of cis-regulatory mechanisms can be
expressed as deep neural networks in which nodes and weights have explicit
physiochemical interpretations. We also demonstrate how such biophysical
networks can be rapidly inferred, using modern deep learning frameworks, from
the data produced by certain types of massively parallel reporter assays
(MPRAs). These results suggest a scalable strategy for using MPRAs to
systematically characterize the biophysical basis of gene regulation in a wide
range of biological contexts. They also highlight gene regulation as a
promising venue for the development of scientifically interpretable approaches
to deep learning.

Deep learning – the use of large multi-layer neural networks in machine
learning applications – is revolutionizing information technology [1]. There
is currently a great deal of interest in applying deep learning techniques to
problems in genomics, especially for understanding gene regulation [2, 3, 4,
5]. These applications remain somewhat controversial, however, due to the
difficulty of mechanistically interpreting neural network models trained on
functional genomics data. Multiple attribution strategies, which seek to
extract meaning post-hoc from neural networks that have rather generic
architectures, have been proposed for addressing this interpretability problem
[6, 7, 8]. However, there remains a substantial gap between the outputs of
such attribution methods and fully mechanistic models of gene regulation.

Here we advocate for a complementary approach: the inference of neural network
models whose architecture reflects explicit biophysical hypotheses for how
cis-regulatory sequences function. This strategy is based on two key
observations. First, two large classes of biophysical models can be formulated
as three-layer neural networks that have a stereotyped form and in which nodes
and weights have explicit physiochemical interpretations. This is true of
thermodynamic models, which rely on a quasi-equilibrium assumption [9, 10, 11,
12, 13, 14], as well as kinetic models, which are more complex but do not make
such assumptions [15, 16, 17]. Second, existing deep learning frameworks such
as TensorFlow [18] are able to rapidly infer such models from the data
produced by certain classes of MPRAs.

The idea of using neural networks to model the biophysics of gene regulation
goes back to at least [19]. To our knowledge, however, the fact that all
thermodynamic models of gene regulation can be represented by a simple three-
layer architecture has not been previously reported. We are not aware of any
prior work that uses neural networks to represent kinetic models or models
involving King-Altman diagrams. There is a growing literature on modeling MPRA
data using neural networks [20, 21, 22, 23, 24]. However, such modeling has
yet to be advocated for MPRAs that dissect the mechanisms of single cis-
regulatory sequences, such as in [25, 26, 27, 28, 29], which is our focus
here. Finally, some recent work has used deep learning frameworks to infer
parametric models of cis-regulation [30, 31]. These studies did not, however,
infer the type of explicit biophysical quantities, such as $\Delta G$ values,
that our approach recovers.

Figure 1: A thermodynamic model of transcriptional regulation. (A)
Transcriptional activation at the E. coli lac promoter is regulated by two
proteins, CRP and $\upsigma^{70}$ RNA polymerase (RNAP). CRP is a
transcriptional activator that up-regulates transcription by stabilizing RNAP-
DNA binding. $\Delta G_{C}$ and $\Delta G_{R}$ respectively denote the Gibbs
free energies of the CRP-DNA and RNAP-DNA interactions, while $\Delta G_{I}$
denotes the Gibbs free energy of interaction between CRP and RNAP. (B) Like
all thermodynamic models of gene regulation, this model consists of a set of
states, each state having an associated Gibbs free energy and activity. The
probability of each state is assumed to follow the Boltzmann distribution. (C)
The corresponding activity predicted by such thermodynamic models is the
state-specific activity averaged together using these Boltzmann probabilities.
(D) This model formulated as a three-layer neural network. First layer nodes
represent interaction energies, second layer nodes represent state
probabilities, and third layer nodes represent transcriptional activity. The
values of weights are indicated; gray lines correspond to zero weights. The
second layer has a softmin activation, while the third has a linear
activation. All thermodynamic models of cis-regulation can be represented
using this general three-layer form.

All thermodynamic models of cis-regulation can be represented as three-layer
neural networks as follows. First one defines a set of molecular complexes, or
“states”, which we index using $s$. Each state has both a Gibbs free energy
$\Delta G_{s}$ and an associated activity $\alpha_{s}$. These energies
determine the probability $P_{s}$ of each state occurring in thermodynamic
equilibrium via the Boltzmann distribution,111To reduce notational burden, all
$\Delta G$ values are assumed to be in thermal units. At $37^{\circ}$C, one
thermal unit is $1\,k_{B}T=0.62~{}\mathrm{kcal/mol}$, where $k_{B}$ is
Boltzmann’s constant and $T$ is temperature.

$P_{s}=\frac{e^{-\Delta G_{s}}}{\sum_{s^{\prime}}e^{-\Delta G_{s^{\prime}}}}.$
(1)

The energy of each state is, in turn, computed using integral combinations of
the individual interaction energies $\Delta G_{j}$ that occur in that state.
We can therefore write $\Delta G_{s}=\sum_{j}\omega_{sj}\Delta G_{j}$, where
$\omega_{sj}$ is the number of times that interaction $j$ occurs in state $s$.
The resulting activity predicted by the model is given by the activities
$\alpha_{s}$ of the individual states averaged over this distribution, i.e.,
$t=\sum_{s}\alpha_{s}P_{s}$.

Fig. 1 illustrates a thermodynamic model for transcriptional activation at the
E. coli lac promoter. This model involves two proteins, CRP and RNAP, as well
as three interaction energies: $\Delta G_{C}$, $\Delta G_{R}$, and $\Delta
G_{I}$. The rate of transcription $t$ is further assumed to be proportional to
the fraction of time that RNAP is bound to DNA (Fig. 1A). This model is
summarized by four different states, two of which lead to transcription and
two of which do not (Fig. 1B). Fig. 1C shows the resulting formula for $t$ in
terms of model parameters. This model is readily formulated as a feed-forward
neural network (Fig. 1D). Indeed, all thermodynamic models of cis-regulation
can be formulated as three-layer neural networks: layer 1 represents molecular
interaction energies, layer 2 (which uses a softmin activation) represents
state probabilities, and layer 3 (using linear activation) represents the
biological activity of interest, which in this case is transcription rate.

Figure 2: Inference of a thermodynamic model from MPRA data. (A) Schematic of
the sort-seq MPRA of [26]. A 75 bp region of the E. coli lac promoter was
mutagenized at 12% per nucleotide. Variant promoters were then used to drive
the expression of GFP. Cells carrying these expression constructs were then
sorted using FACS, and the variant sequences in each bin were sequenced. This
yielded data on about $5\times 10^{4}$ variant promoters across 10 bins. (B)
The neural network from Fig. 1D, but with $\Delta G_{C}$ and $\Delta G_{R}$
expressed as linear functions of the DNA sequence $\vec{x}$, as well as a
dense feed-forward network mapping activity $t$ to bins via a probability
distribution $p(\mathrm{bin}|t)$. Gray lines indicate weights fixed at 0. The
weights linking nodes $P_{3}$ and $P_{4}$ to node $t$ were constrained to have
the same value $t_{sat}$. (C) The parameter values inferred from the MPRA data
of [26]. Shown are the CRP energy matrix $\vec{\theta}_{C}$, the RNAP energy
matrix $\vec{\theta}_{R}$, and the CRP-RNAP interaction energy $\Delta G_{I}$.
Since increasingly negative energy corresponds to stronger binding, they
$y$-axis in the logo plots is inverted. Logos were generated using Logomaker
[32].

We can infer thermodynamic models like these for a cis-regulatory sequence of
interest (the wild-type sequence) from the data produced by an MPRA performed
on an appropriate sequence library [26]. Indeed, a number of MPRAs have been
performed with this explicit purpose in mind [26, 27, 33, 34, 35]. To this
end, such MPRAs are generally performed using libraries that consist of
sequence variants that differ from the wild-type sequence by a small number of
single nucleotide polymorphisms (SNPs). The key modeling assumption that
motivates using libraries of this form is that the assayed sequence variants
will form the same molecular complexes as the wild-type sequence, but with
Gibbs free energies and state activities whose values vary from sequence to
sequence. By contrast, variant libraries that contain insertions, deletions,
or large regions of random DNA (e.g. [20, 21, 23, 24, 31]) are unlikely to
satisfy this modeling assumption.

Fig. 2A summarizes the sort-seq MPRA described in [26]. Lac promoter variants
were used to drive GFP expression in E. coli, cells were sorted into 10 bins
using fluorescence-activated cell sorting, and the variant promoters within
each bin were sequenced. This yielded data comprising about $5\times 10^{4}$
variant lac promoter sequences, each associated with one of 10 bins. The
authors then fit the biophysical model shown in Fig. 1C, but under the
assumption that $\Delta G_{C}=\vec{\theta}_{C}\cdot\vec{x}+\mu_{C}$ and
$\Delta G_{R}=\vec{\theta}_{R}\cdot\vec{x}+\mu_{R}$, where $\vec{x}$ is a one-
hot encoding of promoter DNA sequence.

Here we used TensorFlow to infer the same model formulated as a deep neural
network. Specifically, we augmented the network in Fig. 1D by making $\Delta
G_{C}$ and $\Delta G_{R}$ sequence-dependent as in [26]. To link $t$ to the
MPRA measurements, we introduced a feed-forward network with one hidden layer
and a softmax output layer corresponding to the 10 bins into which cells were
sorted. Model parameters were then fit to the MPRA dataset using stochastic
gradient descent and early stopping. The results agreed well with those
reported in [26]. In particular, the parameters in the energy matrices for CRP
($\vec{\theta}_{C}$) and RNAP ($\vec{\theta}_{R}$) exhibited respective
Pearson correlation coefficients of 0.986 and 0.994 with those reported in
[26]. The protein-protein interaction energy that we found, $\Delta
G_{I}=-2.9$ kcal/mol, was also compatible with the previously reported value
$\Delta G_{I}=-3.3\pm 0.4$ kcal/mol.

A major difference between our results and those of [26] is the ease with
which they were obtained. Training of the network in Fig. 2B consistently took
about 15 minutes on a standard laptop computer. The model fitting procedure in
[26], by contrast, relied on a custom Parallel Tempering Monte Carlo algorithm
that took about a week to run on a multi-node computer cluster (personal
communication), and more recent efforts to train biophysical models on MPRA
data have encountered similar computational bottlenecks [34, 35].

Also of note is the fact that in [26] the authors inferred models using
information maximization. Specifically, the authors fit the parameters of
$t(\vec{x})$ by maximizing the mutual information $I[t;\mathrm{bin}]$ between
model predictions and observed bins. One difficulty with this strategy is the
need to estimate mutual information. Instead, we used maximum likelihood to
infer the parameters of $t(\vec{x})$ as well as the experimental transfer
function (i.e., noise model) $p(\mathrm{bin}|t)$, which was modeled by a dense
feed-forward network with one hidden layer. These two inference methods,
however, are essentially equivalent: in the large data regime, the parameters
of $t$ that maximize $I[t;\mathrm{bin}]$ are the same as the parameters one
obtains when maximizing likelihood over the parameters of both $t$ and
$p(\mathrm{bin}|t)$; see [36, 37, 38].

Figure 3: A kinetic model for transcriptional initiation by E. coli RNAP. (A)
In this model, promoter DNA can participate in three complexes: unbound,
closed, and open [39, 40]. Transitions between these complexes are governed by
four rate constants: $k_{1}$, $k_{-1}$, $k_{2}$, and $k_{3}$. (B) A formula
for the steady-state rate of mRNA production can be obtained using King-Altman
diagrams [41, 42]. (C) This formula can be represented using the three-layer
neural network shown, where layer 1 represents log transition rates, layer 2
(softmax activation) represents normalized King-Altman diagrams, and layer 3
(linear activation) represents promoter activity. Black lines indicate weight
1; gray lines indicate weight 0. Note that the single nonzero weight
connecting layer 2 to layer 3 (orange) is actually the transition rate $k_{3}$
from layer 1. All kinetic models of cis-regulation share this general three-
layer form.

A shortcoming of thermodynamic models is that they ignore non-equilibrium
processes. Kinetic models address this problem by providing a fully non-
equilibrium characterization of steady-state activity. Kinetic models are
specified by listing explicit state-to-state transition rates rather than
Gibbs free energies. Fig. 3A shows a three-state kinetic model of
transcriptional initiation consisting of unbound promoter DNA, an RNAP-DNA
complex in the closed conformation, and the RNAP-DNA complex in the open
conformation [39, 40]. The rate $k_{3}$ going from the open state to the
unbound state represents transcript initiation. The overall transcription rate
in steady state is therefore $k_{3}$ times the occupancy of the open complex.

King-Altman diagrams [41, 42], a technique from mathematical enzymology,
provide a straight-forward way to compute steady-state occupancy in kinetic
models. Specifically, each state’s occupancy is proportional to the sum of
directed spanning trees (a.k.a. King-Altman diagrams) that flow to that state,
where each spanning tree’s value is given by the product of rates comprising
that tree. Fig. 3B illustrates this procedure for the kinetic model in Fig.
3A. Every such kinetic model can be represented by a three-layer neural
network (e.g., Fig. 3C) in which first layer nodes represent log transition
rates, second layer nodes (after a softmax activation) represent normalized
King-Altman diagrams, and third layer nodes represent the activities of
interest.

Here we have shown how both thermodynamic and kinetic models of gene
regulation can be formulated as three-layer deep neural networks in which
nodes and weights have explicit biophysical meaning. This represents a new
strategy for interpretable deep learning in the study of gene regulation, one
complementary to existing post-hoc attribution methods. We have further
demonstrated that a neural-network-based thermodynamic model can be rapidly
inferred from MPRA data using TensorFlow. This was done in the context of a
well-characterized bacterial promoter because previous studies of this system
have established concrete results against which we could compare our inferred
model. The same modeling approach, however, should be readily applicable to a
wide variety of biological systems amenable to MPRAs, including
transcriptional regulation and alternative mRNA splicing in higher eukaryotes.

Code Availability and Acknowledgements: The neural network model shown in Fig.
2C, as well as the scripts used to infer it from the data of [26], are
available at https://github.com/jbkinney/19_mlcb. We thank Anand Murugan,
Yifei Huang, Peter Koo, Alan Moses, and Mahdi Kooshkbaghi for helpful
discussions, as well as and three anonymous referees for providing
constructive criticism. This project was supported by NIH grant 1R35GM133777
and a grant from the CSHL/Northwell Health partnership.

## References

  * [1] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” Nature, vol. 521, pp. 436–444, May 2015.
  * [2] J. Zhou and O. G. Troyanskaya, “Predicting effects of noncoding variants with deep learning–based sequence model,” Nat Methods, vol. 12, pp. 931–934, Aug. 2015.
  * [3] B. Alipanahi, A. Delong, M. T. Weirauch, and B. J. Frey, “Predicting the sequence specificities of dna- and rna-binding proteins by deep learning,” Nat Biotechnol, vol. 33, pp. 831–838, Jul 2015.
  * [4] K. Jaganathan, S. K. Panagiotopoulou, J. F. McRae, S. F. Darbandi, D. Knowles, Y. I. Li, J. A. Kosmicki, J. Arbelaez, W. Cui, G. B. Schwartz, E. D. Chow, E. Kanterakis, H. Gao, A. Kia, S. Batzoglou, S. J. Sanders, and K. K.-H. Farh, “Predicting splicing from primary sequence with deep learning,” Cell, vol. 176, pp. 535–548.e24, Jan. 2019.
  * [5] G. Eraslan, Ž. Avsec, J. Gagneur, and F. J. Theis, “Deep learning: new computational modelling techniques for genomics,” Nat Rev Genet, vol. 20, no. 7, pp. 389–403, 2019.
  * [6] K. Simonyan, A. Vedaldi, and A. Zisserman, “Deep inside convolutional networks: Visualising image classification models and saliency maps,” arXiv preprint arXiv:1312.6034, 2013.
  * [7] A. Shrikumar, P. Greenside, and A. Kundaje, “Learning important features through propagating activation differences,” in Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, pp. 3145–3153, JMLR.org, 2017.
  * [8] A. Chattopadhyay, P. Manupriya, A. Sarkar, and V. N. Balasubramanian, “Neural network attributions: A causal perspective,” arXiv preprint arXiv:1902.02302, 2019.
  * [9] G. Ackers, A. Johnson, and M. Shea, “Quantitative model for gene regulation by lambda phage repressor,” Proc Natl Acad Sci USA, vol. 79, pp. 1129–1133, Feb. 1982.
  * [10] M. A. Shea and G. K. Ackers, “The or control system of bacteriophage lambda. a physical-chemical model for gene regulation,” J Mol Biol, vol. 181, pp. 211–230, Jan. 1985.
  * [11] L. Bintu, N. E. Buchler, H. G. Garcia, U. Gerland, T. Hwa, J. Kondev, T. Kuhlman, and R. Phillips, “Transcriptional regulation by the numbers: applications,” Curr Opin Genet Dev, vol. 15, pp. 125–135, Apr. 2005.
  * [12] L. Bintu, N. E. Buchler, H. G. Garcia, U. Gerland, T. Hwa, J. Kondev, and R. Phillips, “Transcriptional regulation by the numbers: models,” Curr Opin Genet Dev, vol. 15, pp. 116–124, Apr. 2005.
  * [13] E. Segal and J. Widom, “From dna sequence to transcriptional behaviour: a quantitative approach,” Nat Rev Genet, vol. 10, pp. 443–456, July 2009\.
  * [14] M. S. Sherman and B. A. Cohen, “Thermodynamic state ensemble models of cis-regulation,” PLoS Comput Biol, vol. 8, no. 3, p. e1002407, 2012.
  * [15] J. Estrada, F. Wong, A. DePace, and J. Gunawardena, “Information integration and energy expenditure in gene regulation,” Cell, vol. 166, pp. 234–244, June 2016.
  * [16] C. Scholes, A. H. DePace, and A. Sanchez, “Combinatorial gene regulation through kinetic control of the transcription cycle,” Cell Syst, vol. 4, pp. 97–108.e9, Jan. 2017.
  * [17] J. Park, J. Estrada, G. Johnson, B. J. Vincent, C. Ricci-Tam, M. D. Bragdon, Y. Shulgina, A. Cha, Z. Wunderlich, J. Gunawardena, and A. H. DePace, “Dissecting the sharp response of a canonical developmental enhancer reveals multiple sources of cooperativity,” eLife, vol. 8, p. 2787, June 2019.
  * [18] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, et al., “Tensorflow: A system for large-scale machine learning,” in 12th USENIX Symposium on Operating Systems Design and Implementation (OSDI 16), pp. 265–283, 2016.
  * [19] E. Mjolsness, “On cooperative quasi-equilibrium models of transcriptional regulation,” J Bioinform Comput Biol, vol. 5, no. 2b, pp. 467–490, 2007\.
  * [20] A. B. Rosenberg, R. P. Patwardhan, J. Shendure, and G. Seelig, “Learning the sequence determinants of alternative splicing from millions of random sequences,” Cell, vol. 163, pp. 698–711, Oct. 2015.
  * [21] J. T. Cuperus, B. Groves, A. Kuchina, A. B. Rosenberg, N. Jojic, S. Fields, and G. Seelig, “Deep learning of the regulatory grammar of yeast 5’ untranslated regions from 500,000 random sequences,” Genome Res, vol. 27, no. 12, pp. 2015 – 2024, 2017.
  * [22] R. Movva, P. Greenside, G. K. Marinov, S. Nair, A. Shrikumar, and A. Kundaje, “Deciphering regulatory dna sequences and noncoding genetic variants using neural network models of massively parallel reporter assays,” PLoS ONE, vol. 14, no. 6, p. e0218073, 2019.
  * [23] P. J. Sample, B. Wang, D. W. Reid, V. Presnyak, I. J. McFadyen, D. R. Morris, and G. Seelig, “Human 5’ utr design and variant effect prediction from a massively parallel translation assay,” Nat Biotechnol, vol. 37, no. 7, pp. 803–809, 2019.
  * [24] N. Bogard, J. Linder, A. B. Rosenberg, and G. Seelig, “A deep neural network for predicting and engineering alternative polyadenylation,” Cell, vol. 178, no. 1, pp. 91–106.e23, 2019.
  * [25] R. P. Patwardhan, C. Lee, O. Litvin, D. L. Young, D. Pe’er, and J. Shendure, “High-resolution analysis of dna regulatory elements by synthetic saturation mutagenesis,” Nat Biotechnol, vol. 27, no. 12, pp. 1173 – 1175, 2009.
  * [26] J. B. Kinney, A. Murugan, C. G. Callan, and E. C. Cox, “Using deep sequencing to characterize the biophysical mechanism of a transcriptional regulatory sequence,” Proc Natl Acad Sci USA, vol. 107, pp. 9158–9163, May 2010.
  * [27] A. Melnikov, A. Murugan, X. Zhang, T. Tesileanu, L. Wang, P. Rogov, S. Feizi, A. Gnirke, C. G. Callan, J. B. Kinney, M. Kellis, E. S. Lander, and T. S. Mikkelsen, “Systematic dissection and optimization of inducible enhancers in human cells using a massively parallel reporter assay,” Nat Biotechnol, vol. 30, pp. 271–277, Feb. 2012.
  * [28] J. C. Kwasnieski, I. Mogno, C. A. Myers, J. C. Corbo, and B. A. Cohen, “Complex effects of nucleotide variants in a mammalian cis-regulatory element,” Proc Natl Acad Sci USA, vol. 109, pp. 19498–19503, Nov. 2012\.
  * [29] R. P. Patwardhan, J. B. Hiatt, D. M. Witten, M. J. Kim, R. P. Smith, D. May, C. Lee, J. M. Andrie, S.-I. Lee, G. M. Cooper, N. Ahituv, L. A. Pennacchio, and J. Shendure, “Massively parallel functional dissection of mammalian enhancers in vivo,” Nat Biotechnol, vol. 30, no. 3, pp. 265 – 270, 2012\.
  * [30] Y. Liu, K. Barr, and J. Reinitz, “Fully interpretable deep learning model of transcriptional control,” bioRxiv preprint doi:10.1101/655639, May 2019\.
  * [31] C. G. d. Boer, E. D. Vaishnav, R. Sadeh, E. L. Abeyta, N. Friedman, and A. Regev, “Deciphering eukaryotic gene-regulatory logic with 100 million random promoters,” Nat Biotechnol, pp. 1–10, 2019.
  * [32] A. Tareen and J. B. Kinney, “Logomaker: beautiful sequence logos in python,” Bioinformatics, Dec. 2019.  btz921.
  * [33] M. Razo-Mejia, J. Q. Boedicker, D. Jones, A. DeLuna, J. B. Kinney, and R. Phillips, “Comparison of the theoretical and real-world evolutionary potential of a genetic circuit,” Phys Biol, vol. 11, p. 026005, Apr. 2014\.
  * [34] N. M. Belliveau, S. L. Barnes, W. T. Ireland, D. L. Jones, M. J. Sweredoski, A. Moradian, S. Hess, J. B. Kinney, and R. Phillips, “Systematic approach for dissecting the molecular mechanisms of transcriptional regulation in bacteria,” Proc Natl Acad Sci USA, vol. 115, pp. E4796–E4805, May 2018\.
  * [35] S. L. Barnes, N. M. Belliveau, W. T. Ireland, J. B. Kinney, and R. Phillips, “Mapping dna sequence to transcription factor binding energy in vivo,” PLoS Comput Biol, vol. 15, p. e1006226, Feb. 2019.
  * [36] J. B. Kinney, G. Tkačik, and C. G. Callan, “Precise physical models of protein–DNA interaction from high-throughput data,” Proc Natl Acad Sci USA, vol. 104, pp. 501–506, Jan. 2007.
  * [37] J. B. Kinney and G. S. Atwal, “Parametric inference in the large data limit using maximally informative models,” Neural Comput, vol. 26, pp. 637–653, Apr. 2014.
  * [38] G. S. Atwal and J. B. Kinney, “Learning quantitative sequence–function relationships from massively parallel experiments,” J Stat Phys, vol. 162, no. 5, pp. 1203–1243, 2016.
  * [39] W. R. McClure, “Rate-limiting steps in rna chain initiation,” Proc Natl Acad Sci USA, vol. 77, no. 10, pp. 5634 – 5638, 1980.
  * [40] W. R. McClure, “Mechanism and control of transcription initiation in prokaryotes,” Annu Rev Biochem, vol. 54, no. 1, pp. 171 – 204, 1985.
  * [41] E. King and C. Altman, “A schematic method of deriving the rate laws for enzyme-catalyzed reactions,” J Phys Chem, vol. 60, no. 10, pp. 1375–1378, 1956.
  * [42] T. L. Hill, Free Energy Transduction and Biochemical Cycle Kinetics.  New York: Springer-Verlag, 1989.


# Discovering Ancestral Instrumental Variables for Causal Inference from
Observational Data

Debo Cheng and Jiuyong Li and Lin Liu and Kui Yu and Thuc Duy Le and Jixue Liu
D. Cheng, J. Li, L. Liu, T.D. Le and J. Liu are with STEM, University of South
Australia, Mawson Lakes, South Australia, 5095, Australia (e-mail:
{Debo.Cheng,Jiuyong.Li,Lin.Liu,Thuc.Le,Jixue.Liu}@unisa.edu.au). K. Yu is with
the School of Computer Science and Information Engineering, Hefei University
of Technology, Hefei, 230000, China (e-mail: ykui713@gmail.com).
(Corresponding authors: Debo Cheng and Jiuyong Li)

###### Abstract

Instrumental variable (IV) is a powerful approach to inferring the causal
effect of a treatment on an outcome of interest from observational data even
when there exist latent confounders between the treatment and the outcome.
However, existing IV methods require that an IV is selected and justified with
domain knowledge. An invalid IV may lead to biased estimates. Hence,
discovering a valid IV is critical to the applications of IV methods. In this
paper, we study and design a data-driven algorithm to discover valid IVs from
data under mild assumptions. We develop the theory based on partial ancestral
graphs (PAGs) to support the search for a set of candidate Ancestral IVs
(AIVs), and for each possible AIV, the identification of its conditioning set.
Based on the theory, we propose a data-driven algorithm to discover a pair of
IVs from data. The experiments on synthetic and real-world datasets show that
the developed IV discovery algorithm estimates accurate estimates of causal
effects in comparison with the state-of-the-art IV based causal effect
estimators.

###### Index Terms:

Causal Inference, Observational Studies, Latent Confounders, Confounding Bias,
Maximal Ancestral Graph.

##  I Introduction

Afundamental challenge for inferring from observational data the causal effect
of a treatment $W$ (a.k.a. exposure, intervention or action) on an outcome $Y$
of interest is the presence of latent (a.k.a. unobserved or unmeasured)
confounders, variables which affect $W$ and $Y$ simultaneously. Instrumental
variables (IVs) are a powerful tool to address this challenge, primarily used
by statisticians, economists and social scientist [1, 2, 3]. It is possible to
eliminate the confounding bias by leveraging a valid IV [4].

The standard IV approach requires a predefined IV (denoted as $S$) that meets
the following three conditions: (1) $S$ is a cause of $W$, (2) There is no
confounding bias for the effect of $S$ on $Y$ (a.k.a. exogeneity), and (3) the
effect of $S$ on $Y$ is entirely mediated through $W$ (a.k.a. the exclusion
restriction) [5]. For example, $S$ in the causal graph in Fig. 1 (a) is a
standard IV as it meets all the three conditions. The second and third
conditions have to be justified by domain knowledge [6] and hence IV based
methods are mostly classified as empirical methods in literature [7, 8].

Figure 1: (a) $S$ is a standard IV for $W\rightarrow Y$. (b) $S_{1}$ and
$S_{2}$ are valid IVs wrt., $W\rightarrow Y$, conditioning on
$\mathbf{Z}=\\{X_{2},X_{3}\\}$.

Data-driven approach to identifying standard IVs is impractical. Instrumental
inequality [9, 7] has been proposed to test an IV in data, but it is a
necessary condition, not a sufficient condition. Under a set of strong
assumptions on the data distribution, a semi-instrumental variable that can be
tested in data has been introduced by [10]. However, there is not a practical
algorithm rooted from the above concepts. By assuming that at least a half of
the covariates are valid IVs, Kang et al.  [11] proposed an algorithm, sisVIVE
(some invalid, some valid IV estimator), to estimate the causal effect in
data. Hartford et al.  [12] extended the sisVIVE algorithm by employing a deep
learning based IV estimator. The main challenge of the two algorithms is that
their strong assumption is often unsatisfied in many real-world applications.

To relax the last two conditions (i.e., the exogeneity and exclusion
restriction) of a standard IV, a graphical criterion [6, 13] is proposed to
identify an observed variable $S$ as an IV (i.e., the conditional IV (CIV)),
conditioning on a set of observed variables $\mathbf{Z}$ from a given DAG
(directed acyclic graph, which represents the causal relations of all measured
and unmeasured variables). CIV allows a confounding bias between $S$ and $Y$,
and $S$ have multiple causal paths to $Y$. The bias can be adjusted by a
conditioning set. Van der Zander et al.  [14] have revised CIV to Ancestral IV
(AIV) to avoid the situation where void CIVs may be identified based on the
original definition (see Section II-B for details). AIV identification needs a
DAG too.

IV.tetrad [8] is the only existing data-driven CIV method. IV.tetrad based on
CIV and it requires two valid CIVs in the covariate set. The tetrad condition
is used to discover the pair of CIVs in data and it assumes that the
conditioning sets of the pair of CIVs are the same and equal to the set of
remaining covariates (i.e. the original covariate set excluding the pair of
CIVs). This tetrad condition leads to a wrong identification when the set of
remaining covariates contains a collider. For example, $S_{1}$ and $S_{2}$ in
Fig. 1 (b) are a pair of CIVs. The conditioning sets for $S_{1}$ and $S_{2}$
are the same, both equal to $\\{X_{2},X_{3}\\}$, but not as assumed in the
IV.tetrad, i.e. $\\{X_{1},X_{2},X_{3}\\}$. When $X_{1}$ is used in the
conditioning set, path $S_{1}\to X_{1}\leftarrow U_{2}\to Y$ is opened and
$S_{1}$ does not instrumentalise $W$ conditioning on $\\{X_{1},X_{2},X_{3}\\}$
any more. Hence, the tetrad condition does not find the right CIV pairs and
leads to a biased causal effect estimation.

This work improves IV.tetrad in the following ways, which are also our
contributions.

  1. 1.

We generalise the tetrad condition so that each AIV in the pair conditions on
its own conditioning set, and this rectifies the current tetrad condition
which fails to find the right pair when the covariate set contains a collider.

  2. 2.

We develop the theory for identifying the set of candidate AIVs in a reduced
space for efficient search for a pair of AIVs.

  3. 3.

We propose a data-driven algorithm for estimating causal effects from data
with latent variables based on the above developed theorems. Extensive
experiments on synthetic and real-world datasets have shown the effectiveness
of the algorithm.

##  II Background

###  II-A Graph Terminology

A graph $\mathcal{G}\\!\\!=\\!\\!(\mathbf{V},\mathbf{E})$ is composed of a set
of nodes $\mathbf{V}\\!\\!=\\!\\!\\{V_{1},\dots,V_{p}\\}$, representing random
variables, and a set of edges
$\mathbf{E}\\!\subseteq\mathbf{V}\times\mathbf{V}$, representing the relations
between nodes. In this paper, we assume that in a graph $\mathcal{G}$, there
is at most one edge between any two nodes.

Two nodes are adjacent if there exists an edge between them. For an edge
$V_{i}\rightarrow V_{j}$, $V_{i}$ and $V_{j}$ are its head and tail
respectively and $V_{i}$ is known as a parent of $V_{j}$ (and $V_{j}$ is a
child of $V_{i}$). We use $Adj(V)$, $Pa(V)$ and $Ch(V)$ to denote the sets of
all adjacent nodes, parents and children of $V$, respectively. A path $\pi$ is
a sequence of nodes $\left\langle V_{1},\dots,V_{n}\right\rangle$ such that
for $1\leq i\leq n-1$, the pair $(V_{i},V_{i+1})$ is adjacent. A path $\pi$
from $V_{i}$ to $V_{j}$ is a directed or causal path if all edges along it are
directed towards $V_{j}$, and $(V_{i},V_{j})$ are called endpoint nodes, other
nodes are non-endpoint nodes. If there is a directed path $\pi$ from $V_{i}$
to $V_{j}$, $V_{i}$ is a known as an ancestor of $V_{j}$ and $V_{j}$ is a
descendant of $V_{i}$. The sets of ancestors and descendants of a node $V$ are
denoted as $An(V)$ and $De(V)$, respectively.

A DAG (directed acyclic graph) is a direct graph (i.e. a graph containing only
directed edges $\rightarrow$) without directed cycles (i.e. a directed path
whose two endpoints are the same node). A DAG is often used to represent the
data generation mechanism or causal mechanism underlying the data, with all
variables, both observed and unobserved (if any) included in the graph.

Ancestral graphs are used to represent the data generating generation
mechanisms that may involve latent variables, with only observed variables
included in the graphs [15, 16]. An ancestral graph may contain three types of
edges, $\rightarrow$, $\leftrightarrow$ (it is used to represent that a common
cause of two observed variables is a latent variable.) and
$\put(2.0,2.5){\circle{2.5}}\rightarrow$ (the circle tail denotes the
orientation of the edge is uncertain.), and we use
‘$\put(0.2,-2.2){*}\rightarrow$’ to denote any of the three types. $C$ is a
collider on the path $\pi$ if $\pi$ contains a subpath
$\put(0.2,-2.2){*}\rightarrow C\leftarrow\put(-5.0,-2.2){*}$. In an ancestral
graph, a path is a collider path if every non-endpoint node on it is a
collider. A path of length one is a trivial collider path.

In a graph, an almost directed cycle occurs when $V_{i}\leftrightarrow V_{j}$
is in the graph and $V_{j}\in An(V_{i})$. An ancestral graph is a graph that
does not contain directed cycles or almost directed cycles [15]. In an
ancestral graph, a path from $V_{i}$ to $V_{j}$ is a possibly directed or
causal path if there is not an arrowhead pointing in the direction of $V_{i}$.
In this case, $V_{i}$ is a possible ancestor of $V_{j}$ and $V_{j}$ is a
possible descendant of $V_{i}$. The sets of possible ancestors and descendants
of $V$ are denoted as $PossAn(V)$ and $PossDe(V)$, respectively.

In graphical causal modelling, the assumptions of Markov property and
faithfulness are often involved to discuss the relationship between the causal
graph and the data distribution.

######  Definition 1 (Markov property [6]).

Given a DAG $\mathcal{G}=(\mathbf{V},\mathbf{E})$ and the joint probability
distribution of $\mathbf{V}$ $(prob(\mathbf{V}))$, $\mathcal{G}$ satisfies the
Markov property if for $\forall V_{i}\in\mathbf{V}$, $V_{i}$ is
probabilistically independent of all of its non-descendants, given
$Pa(V_{i})$.

######  Definition 2 (Faithfulness [17]).

A DAG $\mathcal{G}=(\mathbf{V},\mathbf{E})$ is faithful to a joint
distribution $prob(\mathbf{V})$ over the set of variables $\mathbf{V}$ if and
only if every independence present in $prob(\mathbf{V})$ is entailed by
$\mathcal{G}$ and satisfies the Markov property. A joint distribution
$prob(\mathbf{V})$ over the set of variables $\mathbf{V}$ is faithful to the
DAG $\mathcal{G}$ if and only if the DAG $\mathcal{G}$ is faithful to the
joint distribution $prob(\mathbf{V})$.

A DAG $\mathcal{G}=(\mathbf{V},\mathbf{E})$ satisfies both assumptions, the
probability distribution $prob(\mathbf{V})$ can be factorised as:
$prob(\mathbf{V})=\prod_{i}^{p}prob(V_{i}|Pa(V_{i}))$. Thus, together Markov
property and faithfulness establish a close relation between the causal graph
and the data distribution.

######  Definition 3 (Causal sufficiency [17]).

In a data, for every pair of observed variables $(V_{i},V_{j})$ in
$\mathbf{V}$, all their common causes are also in $\mathbf{V}$.

In a DAG, d-separation is a well-known graphical criterion that is used to
read off the identification of conditional independence between variables
entailed in the DAG when the Markov property, faithfulness and causal
sufficiency are satisfied [6, 17].

######  Definition 4 (d-separation [6]).

A path $\pi$ in a DAG $\mathcal{G}=(\mathbf{V},\mathbf{E})$ is said to be
d-separated (or blocked) by a set of nodes $\mathbf{Z}$ if and only if (i)
$\pi$ contains a chain $V_{i}\rightarrow V_{k}\rightarrow V_{j}$ or a fork
$V_{i}\leftarrow V_{k}\rightarrow V_{j}$ such that the middle node $V_{k}$ is
in $\mathbf{Z}$, or (ii) $\pi$ contains a collider $V_{k}$ such that $V_{k}$
is not in $\mathbf{Z}$ and no descendant of $V_{k}$ is in $\mathbf{Z}$. A set
$\mathbf{Z}$ is said to d-separate $V_{i}$ from $V_{j}$
($V_{i}\Vbar_{d}V_{j}\mid\mathbf{Z}$) if and only if $\mathbf{Z}$ blocks every
path between $V_{i}$ to $V_{j}$. Otherwise they are said to be d-connected by
$\mathbf{Z}$, denoted as $V_{i}\nVbar_{d}V_{j}\mid\mathbf{Z}$.

######  Property 1.

Two observed variables $V_{i}$ and $V_{j}$ are d-separated given a
conditioning set $\mathbf{Z}$ in a DAG if and only if $V_{i}$ and $V_{j}$ are
conditionally independence given $\mathbf{Z}$ in data [17]. if $V_{i}$ and
$V_{j}$ are d-connected, $V_{i}$ and $V_{j}$ are conditionally dependent.

However, a system may involve the latent variables (the latent variable is an
unmeasured common cause of two nodes) in most situations since there is not a
close world. Ancestral graphs are proposed to represent the system that may
involve latent variables [15, 16]. In our work, we utilise the Maximal
ancestral graph and introduce it as follows.

######  Definition 5 (Maximal ancestral graph (MAG) [15]).

An ancestral graph $\mathcal{M}=(\mathbf{V},\mathbf{E})$ is a MAG when every
pair of non-adjacent nodes $V_{i}$ and $V_{j}$ in $\mathcal{M}$ is m-separated
by a set $\mathbf{Z}\subseteq\mathbf{V}\backslash\\{V_{i},V_{j}\\}$.

It is worth noting that a DAG satisfies both conditions of a MAG, so a DAG is
also a MAG without bi-directed edges [16]. An important concept in a DAG is
d-separation, which captures the conditional independence relationships
between variables based on Markov property [6]. A natural extension of the
d-separation to an ancestral graph is m-separation [15].

######  Definition 6 (m-separation [15]).

In an ancestral graph $\mathcal{M}=(\mathbf{V},\mathbf{E})$, a path $\pi$
between $V_{i}$ and $V_{j}$ is said to be m-separated by a set of nodes
$\mathbf{Z}\subseteq\mathbf{V}\setminus\\{V_{i},V_{j}\\}$ (possibly
$\emptyset$) if $\pi$ contains a subpath $\langle V_{l},V_{k},V_{s}\rangle$
such that the middle node $V_{k}$ is a non-collider on $\pi$ and
$V_{k}\in\mathbf{Z}$; or $\pi$ contains $V_{l}\put(0.2,-2.2){*}\rightarrow
V_{k}\leftarrow\put(-5.0,-2.2){*}V_{s}$ such that $V_{k}\notin\mathbf{Z}$ and
no descendant of $V_{k}$ is in $\mathbf{Z}$. Two nodes $V_{i}$ and $V_{j}$ are
said to be m-separated by $\mathbf{Z}$ in $\mathcal{M}$, denoted as
$V_{i}\Vbar_{m}V_{j}|\mathbf{Z}$ if every path between $V_{i}$ and $V_{j}$ are
m-separated by $\mathbf{Z}$; otherwise they are said to be m-connected by
$\mathbf{Z}$, denoted as $V_{i}\nVbar_{m}V_{j}|\mathbf{Z}$.

where $\Vbar_{m}$ denotes m-separation and $\nVbar_{m}$ denotes m-connecting.
In a DAG, m-separation reduces to d-separation. The Markov property of
ancestral graph is captured by m-separation.

If two MAGs represent the same set of m-separations, they are called _Markov
equivalent_ , and formally, we have the following definition.

######  Definition 7 (Markov equivalent MAGs [18]).

Two MAGs $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ with the same nodes are said
to be _Markov equivalent_ , denoted $\mathcal{M}_{1}\sim\mathcal{M}_{2}$, if
for all triple nodes $X$, $Y$, $Z$, $X$ and $Y$ are m-separated by $Z$ in
$\mathcal{M}_{1}$ if and only if $X$ and $Y$ are m-separated by $Z$ in
$\mathcal{M}_{2}$.

A set of Markov equivalent MAGs can be encoded uniquely by a _partial
ancestral graph_ (PAG) [16].

######  Definition 8 (PAG [16]).

Let $[\mathcal{M}]$ be the Markov equivalence class of a MAG $\mathcal{M}$.
The PAG $\mathcal{P}$ for $[\mathcal{M}]$ is a partial mixed graph such that
(i). $\mathcal{P}$ has the same adjacent relations among nodes as
$\mathcal{M}$ does; (ii). For an edge, its mark of arrowhead or mark of the
tail is in $\mathcal{P}$ if and only if the same mark of arrowhead or the same
mark of the tail is shared by all MAGs in $[\mathcal{M}]$.

######  Definition 9 (Visibility [16]).

Given a MAG $\mathcal{M}=(\mathbf{V},\mathbf{E})$, a directed edge
$V_{i}\rightarrow V_{j}$ is visible if there is a node $V_{k}\notin
Adj(V_{j})$, such that either there is an edge between $V_{k}$ and $V_{i}$
that is into $V_{i}$, or there is a collider path between $V_{k}$ and $V_{i}$
that is into $V_{i}$ and every node on the path is a parent of $V_{j}$.
Otherwise, $V_{i}\rightarrow V_{j}$ is invisible.

The visible edge is a critical concept in a MAG [16, 19]. A DAG over measured
and unmeasured variables can be mapped to a MAG with measured variables. From
a DAG over $\mathbf{X}\cup\mathbf{U}$ where $\mathbf{X}$ is a set of measured
variables and $\mathbf{U}$ is a set of unmeasured variables, following the
construction rule specified in [18], one can construct a MAG with nodes
$\mathbf{X}$ such that all the conditional independence relationships among
the measured variables entailed by the DAG are entailed by the MAG and vice
versa, and the ancestral relationships in the DAG are maintained in the MAG.

###  II-B Instrumental Variable (IV)

Let $W$ be the treatment variable, $Y$ the outcome, and $\mathbf{X}$ be the
set of all other variables. As in the literature, we consider that
$\mathbf{X}$ contains pretreatment variables only, i.e., for any
$X\in\mathbf{X}$, $X\notin(De(W)\cup De(Y))$ [20, 8, 21]. When estimating the
average causal effect of $W$ on $Y$, denoted as $\beta_{wy}$ from data, we
follow the convention in causal inference literature, that is, the data
distribution is said to be compatible with the underlying causal DAG
$\mathcal{G}$, i.e., the assumptions of Markov property and faithfulness are
satisfied.

The goal of this work is to quantify the average causal effect of $W$ on $Y$,
i.e., $\beta_{wy}$, even when there exist unmeasured variables between $W$ and
$Y$, based on observational data, by extending the existing IV techniques.

In this section, we introduce the background information related to IVs.

######  Definition 10 (Standard IV [22, 23]).

A variable $S$ is said to be an IV wrt., $W\rightarrow Y$, if (i) $S$ has a
causal effect on $W$, (ii) $S$ affects $Y$ only through $W$ (i.e., $S$ has no
direct effect on $Y$), and (iii) $S$ does not share common causes with $Y$.

The last two conditions of a standard IV are untestable and strict. In
practice, $S$ may have other causal paths to $Y$, and $S$ is often confounded
with $Y$ by other measured variables. The concept of conditional IV (CIV) in
DAG is proposed to relax the conditions of a standard IV.

######  Definition 11 (Conditional IV (CIV) in DAG [6, 13]).

Given a DAG $\mathcal{G}=(\mathbf{X}\cup\mathbf{U}\cup\\{W,Y\\},\mathbf{E})$
where $\mathbf{X}$ and $\mathbf{U}$ are measured and unmeasured variables
respectively. A variable $S\in\mathbf{X}$ is said to be a CIV wrt.,
$W\rightarrow Y$, if there exists a set of measured variables
$\mathbf{Z}\subseteq\mathbf{X}$ such that (i) $S\nVbar_{d}W\mid\mathbf{Z}$,
(ii) $S\Vbar_{d}Y\mid\mathbf{Z}$ in $\mathcal{G}_{\underline{W}}$, and (iii)
$\forall Z\in\mathbf{Z}$, $Z\notin De(Y)$ where $\mathcal{G}_{\underline{W}}$
is obtained by removing $W\rightarrow Y$ from $\mathcal{G}$.

For a CIV $S$ as defined in Definition 11, $\mathbf{Z}$ is known to
instrumentalise $S$ in the given DAG. However, a variable may be a CIV when it
has zero causal effect on $W$, and this might result in a misleading
conclusion. To mitigate this issue, Ancestral IV (AIV) in DAG is proposed.

######  Definition 12 (Ancestral IV (AIV) in DAG [14]).

Given a DAG $\mathcal{G}=(\mathbf{X}\cup\mathbf{U}\cup\\{W,Y\\},\mathbf{E})$
where $\mathbf{X}$ and $\mathbf{U}$ are measured and unmeasured variables
respectively. A variable $S\in\mathbf{X}$ is said to be an AIV wrt.,
$W\rightarrow Y$, if there exists a set of measured variables
$\mathbf{Z}\subseteq\mathbf{X}\setminus\\{S\\}$ such that (i)
$S\nVbar_{d}W\mid\mathbf{Z}$, (ii) $S\Vbar_{d}Y\mid\mathbf{Z}$ in
$\mathcal{G}_{\underline{W}}$, and (iii) $\mathbf{Z}$ consists of $An(Y)$ or
$An(S)$ or both and $\forall Z\in\mathbf{Z}$, $Z\notin De(Y)$.

An AIV in DAG is a CIV in DAG, but a CIV may not be an AIV. AIV is a
restricted version of CIV [14]. However, the applications of standard IV, CIV
and AIV are established in a complete causal DAG $\mathcal{G}$, which greatly
limits their capacity in real-world applications.

Recently, Cheng et al. proposed the concept of AIV in MAG and the theorem for
identifying a conditioning set $\mathbf{Z}$ that instrumentalises a given AIV
in PAG [24]. However, the work by Cheng et al. assume an AIV in MAG has been
given, and the focus is on finding the conditioning set for the AIV in MAG
from data. Our work in this paper aims to find an AIV in MAG from data, which
makes it possible for complete data-driven search for AIVs.

######  Proposition 1 (AIV in MAG [24]).

Given a DAG
$\mathcal{G}=(\mathbf{X}\cup\mathbf{U}\cup\\{W,Y\\},\mathbf{E}^{\prime})$ with
the edges $W\rightarrow Y$ and $W\leftarrow U\rightarrow Y$ in
$\mathbf{E}^{\prime}$, and $U\in\mathbf{U}$, and the MAG
$\mathcal{M}=(\mathbf{X}\cup\\{W,Y\\},\mathbf{E})$ is mapped from
$\mathcal{G}$ based on the construction rules [18]. Then if $S$ is an AIV
conditioning on a set of measured variables
$\mathbf{Z}\subseteq\mathbf{X}\setminus\\{S\\}$ in $\mathcal{G}$, $S$ is an
AIV conditioning on a set of measured variables
$\mathbf{Z}\subseteq\mathbf{X}\setminus\\{S\\}$ in $\mathcal{M}$.

######  Theorem 1 (Conditioning set of a given AIV in PAG [24]).

Given a DAG
$\mathcal{G}=(\mathbf{X}\cup\mathbf{U}\cup\\{W,Y\\},\mathbf{E}^{\prime})$ with
the edges $W\rightarrow Y$ and $W\leftarrow U\rightarrow Y$ in
$\mathbf{E}^{\prime}$, and $U\in\mathbf{U}$, and let
$\mathcal{M}=(\mathbf{X}\cup\\{W,Y\\},\mathbf{E})$ be the MAG mapped from
$\mathcal{G}$. From data, the mapped MAG $\mathcal{M}$ is represented by a PAG
$\mathcal{P}=(\mathbf{X}\cup\\{W,Y\\},\mathbf{E}^{\prime\prime})$. For a given
ancestral IV $S$ which is a cause or spouse of $W$, the set $PossAn(S\cup
Y)\setminus\\{W,S\\}$ in the learned $\mathcal{P}$ is a set that
instumentalises $S$ in the DAG $\mathcal{G}$.

It is worth noting that Theorem 1 and the work in [24] are to discover
conditioning set for a given AIV, rather than discovering AIVs and the
conditioning set simultaneously. Hence, discovering an AIV and its
conditioning set simultaneously from data remains unresolved, and it is the
problem to be addressed in this work.

##  III Discovering AIVs based on Graphical Causal Modeling

In this section, we first introduce the generalised tetrad condition. Next, we
propose a set of candidates AIVs in MAG. Then, we propose a theorem to
guarantee that the generalised tetrad condition can be used to discover valid
AIVs from data if there exists a pair of AIVs. Finally, we develop a practical
data-driven algorithm for estimating $\beta_{wy}$ from data.

###  III-A The Generalised tetrad Condition

Let $S_{i}$ and $S_{j}$ be a pair of CIVs given the conditioning set
$\mathbf{X}\setminus\\{S_{i},S_{j}\\}$. Let $\sigma_{s_{i}*y*\mathbf{z}}$
$(\sigma_{s_{j}*y*\mathbf{z}})$ denote the partial covariance of $S_{i}$
($S_{j}$) and $Y$ given $\mathbf{Z}$, and $\sigma_{s_{i}*w*\mathbf{z}}$
$(\sigma_{s_{j}*w*\mathbf{z}})$ denote the partial covariance of $S_{i}$
($S_{j}$) and $W$ given $\mathbf{Z}$. Then, we have
$\beta_{wy}=\sigma_{s_{i}*y*\mathbf{z}}/\sigma_{s_{i}*w*\mathbf{z}}=\sigma_{s_{j}*y*\mathbf{z}}/\sigma_{s_{j}*w*\mathbf{z}}$,
which gives us the following tetrad condition [8]:

$\sigma_{s_{i}*y*\mathbf{z}}\sigma_{s_{j}*w*\mathbf{z}}-\sigma_{s_{i}*w*\mathbf{z}}\sigma_{s_{j}*y*\mathbf{z}}=0$
(1)

The tetrad condition can be tested from data directly. It is a necessary
condition for discovering valid CIVs, which means a pair of variables that are
not valid CIVs can also satisfy the tetrad condition.

We consider a more general setting, where a pair of AIVs $S_{i}$ and $S_{j}$
have different conditioning sets
$\mathbf{Z}_{i}\subseteq\mathbf{X}\setminus\\{S_{i}\\}$ and
$\mathbf{Z}_{j}\subseteq\mathbf{X}\setminus\\{S_{j}\\}$ respectively, and
$\mathbf{Z}_{i}$ and $\mathbf{Z}_{j}$ do not need to be equal. Let
$\sigma_{s_{i}*y*\mathbf{z}_{i}}$ $(\sigma_{s_{j}*y*\mathbf{z}_{j}})$ denote
the partial covariance of $S_{i}$ ($S_{j}$) and $Y$ given $\mathbf{Z}_{i}$
($\mathbf{Z}_{j}$), and $\sigma_{s_{i}*w*\mathbf{z}_{i}}$
($\sigma_{s_{j}*w*\mathbf{z}_{j}}$) denote the partial covariance of $S_{i}$
($S_{j}$) and $W$ given $\mathbf{Z}_{i}$ ($\mathbf{Z}_{j}$). Then we have
$\beta_{wy}=\sigma_{s_{i}*y*\mathbf{z}_{i}}/\sigma_{s_{i}*w*\mathbf{z}_{i}}=\sigma_{s_{j}*y*\mathbf{z}_{j}}/\sigma_{s_{j}*w*\mathbf{z}_{j}}$,
which gives us the following generalised tetrad condition:

$\sigma_{s_{i}*y*\mathbf{z}_{i}}\sigma_{s_{j}*w*\mathbf{z}_{j}}-\sigma_{s_{i}*w*\mathbf{z}_{i}}\sigma_{s_{j}*y*\mathbf{z}_{j}}=0$
(2)

In the following, we will show that the generalised tetrad condition can be
used for finding a pair of AIVs in data if there exists a pair of AIVs. In
comparison with the tetrad condition used in IV.tetrad [8], the search space
of the generalised tetrad condition is larger since
$S_{i},S_{j},\mathbf{Z}_{i}$ and $\mathbf{Z}_{j}$ all vary. In the next
section, we will develop a lemma to reduce the search space.

###  III-B The Theory for Discovering AIVs in MAG

We aim to develop a practical solution for discovering AIVs directly from data
by leveraging the property of a MAG.

We first categorise AIVs into direct AIVs and indirect AIVs. When $S$ is an
AIV and it is an adjacent node of the treatment in the given DAG, its
ancestral or adjacent nodes may be AIVs. We call $S$ a direct AIV and the AIVs
which are ancestral or adjacent nodes of $S$ indirect AIVs if they satisfy
Definition 12. We consider direct AIVs since, in practice, indirect AIVs are
rare. Importantly, direct AIVs have a property to support the data-driven
search for their conditioning sets. In the following discussions, all AIVs are
direct AIVs. We have the following conclusion for discovering the direct AIVs
in a MAG $\mathcal{M}$.

######  Lemma 1 (A direct AIV in MAG).

Given a DAG
$\mathcal{G}=(\mathbf{X}\cup\mathbf{U}\cup\\{W,Y\\},\mathbf{E}^{\prime})$ with
the edges $W\rightarrow Y$ and $W\leftarrow U\rightarrow Y$ in
$\mathbf{E}^{\prime}$, and $U\in\mathbf{U}$, and let
$\mathcal{M}=(\mathbf{X}\cup\\{W,Y\\},\mathbf{E})$ be the MAG mapped from
$\mathcal{G}$. If $S$ is a direct AIV in the DAG $\mathcal{G}$, then $S\in
Adj(Y)\setminus\\{W\\}$ in MAG $\mathcal{M}$.

###### Proof.

Firstly, the edges $W\rightarrow Y$ and $W\leftarrow U\rightarrow Y$ in the
given DAG are represent by an invisible edge $W\rightarrow Y$ in the mapped
MAG $\mathcal{M}$ [18]. For an $S\in\mathbf{X}$ to be an eligible direct AIV
in the DAG $\mathcal{G}$, there are only two cases in the mapped MAG
$\mathcal{M}$. The first case is that $S$ has an edge
$S\put(0.2,-2.2){*}\rightarrow W$, then $S$ must have an edge into $Y$, i.e.,
$S\in Adj(Y)\setminus\\{W\\}$, since otherwise $W\rightarrow Y$ in
$\mathcal{M}$ is visible, which contradicts the invisible edge $W\rightarrow
Y$ in $\mathcal{M}$. The second case is that $S$ has a collider path into $W$
and every collider on the path is in $Pa(Y)$, i.e., $S\nVbar W\mid
Pa(Y)\setminus\\{W\\}$, then $S$ must have an edge into $Y$, i.e., $S\in
Adj(Y)\setminus\\{W\\}$, since otherwise $W\rightarrow Y$ in $\mathcal{M}$ is
visible, which contradicts the invisible edge $W\rightarrow Y$ in
$\mathcal{M}$. Therefore, all direct AIVs in the DAG $\mathcal{G}$ are
included in the set $\mathbf{S}=Adj(Y)\setminus\\{Y\\}$ of the mapped MAG
$\mathcal{M}$. ∎

Lemma 1 provides a set of candidate direct AIVs $Adj(Y)\setminus\\{W\\}$ and
reduces the search space of a direct AIV, i.e., the search space of a direct
AIV is reduced from $\mathbf{O}({|\mathbf{X}|})$ to
$\mathbf{O}({|Adj(Y)\setminus\\{W\\}|})$ where
$|Adj(Y)\setminus\\{W\\}|\ll|\mathbf{X}|$.

Next, we will develop a theorem to show that the generalised tetrad condition
in Eq.(2) can be used to discover a pair of direct AIVs directly from data
with latent variables.

######  Theorem 2.

Given a DAG
$\mathcal{G}=(\mathbf{X}\cup\mathbf{U}\cup\\{W,Y\\},\mathbf{E}^{\prime})$ with
the edges $W\rightarrow Y$ and $W\leftarrow U\rightarrow Y$ in
$\mathbf{E}^{\prime}$, and $U\in\mathbf{U}$, and let
$\mathcal{M}=(\mathbf{X}\cup\\{W,Y\\},\mathbf{E})$ be the MAG mapped from
$\mathcal{G}$. Let
$\mathcal{P}=(\mathbf{X}\cup\\{W,Y\\},\mathbf{E}^{\prime\prime})$ be the PAG
which encodes the set of MAGs Markov equivalent to $\mathcal{M}$. If there
exists a pair of direct AIVs $\\{S_{i},S_{j}\\}\subseteq\mathbf{X}$ in the DAG
$\mathcal{G}$, then $\\{S_{i},S_{j}\\}$ must be in $Adj(Y)\setminus\\{W\\}$ in
the PAG $\mathcal{P}$. Moreover, the two sets,
$\mathbf{Z}_{i}=possAn(S_{i}\cup Y)\setminus\\{W,S_{i}\\}$ and
$\mathbf{Z}_{j}=possAn(S_{j}\cup Y)\setminus\\{W,S_{j}\\}$ in the PAG
$\mathcal{P}$ instrumentalise $S_{i}$ and $S_{j}$ in the DAG $\mathcal{G}$,
respectively. Hence, Eq.(2) holds for $S_{i}$ and $S_{j}$ and their
conditioning sets $\mathbf{Z}_{i}$ and $\mathbf{Z}_{j}$.

###### Proof.

According to Lemma 1, $Adj(Y)\setminus\\{W\\}$ in the mapped MAG $\mathcal{M}$
is the set of candidate direct AIVs in the DAG $\mathcal{G}$. Hence, the set
$Adj(Y)\setminus\\{W\\}$ in the PAG $\mathcal{P}$ must be the set of candidate
direct AIVs because the mapped MAG $\mathcal{M}$ is encoded in the PAG
$\mathcal{P}$. Thus, if $S\in\mathbf{X}$ is a direct AIV in the DAG
$\mathcal{G}$, then $S\in Adj(Y)\setminus\\{W\\}$ in the PAG $\mathcal{P}$
i.e., $\\{S_{i},S_{j}\\}\subseteq\mathbf{S}$ holds. According to Theorem 1,
$\mathbf{Z}_{i}=possAn(S_{i}\cup Y)\setminus\\{W,S_{i}\\}$ and
$\mathbf{Z}_{j}=possAn(S_{j}\cup Y)\setminus\\{W,S_{j}\\}$ in the PAG
$\mathcal{P}$ instrumentalise $S_{i}$ and $S_{j}$ in the DAG $\mathcal{G}$,
respectively. Thus, we have
$\beta_{wy}=\sigma_{s_{i}*y*\mathbf{z}_{i}}/\sigma_{s_{i}*w*\mathbf{z}_{i}}=\sigma_{s_{j}*y*\mathbf{z}_{j}}/\sigma_{s_{j}*w*\mathbf{z}_{j}}$.
Therefore,
$\sigma_{s_{i}*y*\mathbf{z}_{i}}\sigma_{s_{j}*w*\mathbf{z}_{j}}-\sigma_{s_{i}*w*\mathbf{z}_{i}}\sigma_{s_{j}*y*\mathbf{z}_{j}}=0$,
i.e., Eq.(2) holds. ∎

Theorem 2 supports a data-driven algorithm to discover a pair of direct AIVs
$\\{S_{i},S_{j}\\}$ and their corresponding conditioning sets $\mathbf{Z}_{i}$
and $\mathbf{Z}_{j}$ by utilising the generalised tetrad condition. In the
next section, based on the theorem, we will propose a practical algorithm for
estimating $\beta_{wy}$ from data with latent variables.

Note that a significant number of direct AIVs are in both
$Adj(W)\setminus\\{Y\\}$ and $Adj(Y)\setminus\\{W\\}$ in a MAG. Sometimes,
they may be missed from $Adj(Y)\setminus\\{W\\}$ due to the false discoveries
of the structure learning algorithm used [25, 26]. In the corresponding DAG,
the direct AIVs are closer to $W$ than $Y$. To avoid the random fluctuations
without sacrificing much efficiency, in our developed practical algorithm, we
extend the search space of Lemma 1 to $Adj(W\cup Y)\setminus\\{W,Y\\}$. This
only adds minor additional costs to the search process.

###  III-C A Practical Algorithm for Estimating $\beta_{wy}$

We develop a practical data-driven algorithm, AIV.GT (Ancestral IV based on
Generalised Tetrad condition), for estimating $\beta_{wy}$ from data with
latent variables. The pseudocode of AIV.GT is listed in Algorithm 1.

Algorithm 1 AIVs based on the Generalised Tetrad condition (AIV.GT)

1: Input: The set of pretreatment variables $\mathbf{X}$, the treatment $W$,
outcome $Y$ and the dataset $\mathcal{D}$; $\alpha$ =0.05

2: Output: $\hat{\beta}_{wy}$, the causal effect of $W$ on $Y$, or _NA_ , i.e.
lacking knowledge

3: Recover a PAG $\mathcal{P}$ from $\mathcal{D}$ by using the rfci algorithm

4: Obtain $\mathbf{S}=Adj(W\cup Y)\setminus\\{W,Y\\}$ from $\mathcal{P}$

5: if $\left|\mathbf{S}\right|\leqslant 1$ then

6: return _NA_

7: else

8: for each $S_{i}\in\mathbf{S}$ do

9: $\mathbf{Z}_{i}\leftarrow PossAn(S_{i}\cup Y)\setminus\\{S_{i},W,Y\\}$

10: $\hat{\beta}_{i}\leftarrow TSLS(W,Y,S_{i},\mathbf{Z}_{i},\mathcal{D})$

11: end for

12: Initialise $\mathbf{Q}=\emptyset$

13: for each pair $(S_{i},S_{j})\in\mathbf{S}$ do

14: if
$Test.tetrad(W,Y,S_{i},S_{j},\mathbf{Z}_{i},\mathbf{Z}_{j},\mathcal{D},\alpha)$
then

15:
$\epsilon_{ij}=\left|\sigma_{s_{i}*y*\mathbf{z}_{i}}\sigma_{s_{j}*w*\mathbf{z}_{j}}-\sigma_{s_{i}*w*\mathbf{z}_{i}}\sigma_{s_{j}*y*\mathbf{z}_{j}}\right|$

16: $\lambda_{ij}=\left|\epsilon_{ij}-\delta_{ij}\right|$ where
$\delta_{ij}=\left|\hat{\beta}_{i}-\hat{\beta}_{j}\right|$

17: $\mathbf{Q}\leftarrow\mathbf{Q}\cup\lambda_{ij}$

18: end if

19: end for

20: if $\left|\mathbf{Q}\right|$ = $\emptyset$ then

21: return _NA_

22: else

23: return $\hat{\beta}_{wy}=mean(\hat{\beta}_{i},\hat{\beta}_{j})$ where the
consistent score $\lambda_{ij}$ is the smallest in $\mathbf{Q}$

24: end if

25: end if

AIV.GT aims to search for the pair of AIVs from data directly without domain
knowledge. The generalised tetrad condition in Eq.(2) held by a pair of AIVs
and their conditioning sets as described in Theorem 2 if there is a pair of
IVs $\\{S_{i},S_{j}\\}$ in data.

To obtain a reliable result, we propose a consistency score to assess which
paired variables are the most likely AIVs based on the generalised tetrad
condition. Let
$\epsilon_{ij}=\left|\sigma_{s_{i}*y*\mathbf{z}_{i}}\sigma_{s_{j}*w*\mathbf{z}_{j}}-\sigma_{s_{i}*w*\mathbf{z}_{i}}\sigma_{s_{j}*y*\mathbf{z}_{j}}\right|$,
and $\delta_{ij}=\left|\hat{\beta}_{i}-\hat{\beta}_{j}\right|$ where
$\hat{\beta}_{i}$ and $\hat{\beta}_{j}$ are the estimated causal effects of
$W$ on $Y$ by using $S_{i}$ and $S_{j}$ as an instrument, respectively. The
consistency score is defined as
$\lambda_{ij}=\left|\epsilon_{ij}-\delta_{ij}\right|$.

The justification of the consistency score is that $\epsilon_{ij}$ is expected
to be close to 0, and the same with $\delta_{ij}$ 0 if the variables $S_{i}$
and $S_{j}$ are AIVs. Theoretically, the pair of variables with either the
smallest $\epsilon_{ij}$ or $\delta_{ij}$ are most likely to be the pair of
IVs, but in practical cases, a pair of variables passing the generalised
tetrad condition test (i.e., a small enough $\epsilon_{ij}$), may have a large
$\delta_{ij}$, or vice versa, since the pair of variables are not IVs. So we
use their difference $\lambda_{ij}$ to avoid such cases because $\lambda_{ij}$
must be smaller than both $\epsilon_{ij}$ and $\delta_{ij}$. Under the
assumption that there exists at least a pair of IVs, if the consistency score
of a pair of variables is the smallest, then the pair are most likely to be
AIVs. The paired variables with the minimal consistency score is returned as
the result of AIV.GT.

The AIV.GT algorithm is divided into two parts. The first part (Lines 3 to 11)
is to obtain all candidate AIV pairs and the possible causal effects of $W$ on
$Y$ estimated using these pairs. Line 3 aims to learn a PAG $\mathcal{P}$ from
data by using a causal structure learning algorithm [27, 26]. We use rfci
(really fast causal inference) [28] in AIV.GT. Line 4 aims to get the set of
candidate AIVs from the learned PAG $\mathcal{P}$. Line 5 tests the size of
$\mathbf{S}$ and if $\left|\mathbf{S}\right|\leqslant 1$, then AIV.GT returns
NA due to lack of knowledge. Lines 8 to 11 are to estimate the causal effect
$\hat{\beta}_{wy}$ using each candidate AIV. Line 9 is to find the
conditioning set $\mathbf{Z}_{i}$ for a candidate AIV $S_{i}$ based on Theorem
1. Line 10, the function $TSLS()$ is the estimator of two-stage least squares
(TSLS) by using $S_{i}$ as an IV and conditioning on $\mathbf{Z}_{i}$ for
calculating $\hat{\beta}_{i}$.

The second part of AIV.GT is to discover the pair of AIVs. Line 12 is to
initialise the set of consistency scores $\mathbf{Q}$. Lines 13 to 19 are to
check the validity of each pair of candidate AIVs based on Theorem 2. If the
generalised tetrad condition holds on a pair of candidate AIVs, then calculate
their consistency score. Line 14, the function $Test.tetrad()$ is implemented
by using the Wishart test wrt., the generalised tetrad condition [29, 17].
$Test.tetrad()$ returns TRUE if and only if the set of candidate pair
variables returns a p-value greater than the significant level $\alpha$
($\alpha$ =0.05 in this work). Lines 15 to 17 are to obtain the consistency
score of each paired AIVs satisfying the generalised tetrad condition. In
Lines 20 and 21, if $\mathbf{Q}$ is an empty set, then no pair of variables
has passed the tetrad condition test and the algorithm returns NA. In Lines 22
to 24, AIV.GT returns the mean causal effect of the pair of variables with the
smallest $\lambda_{ij}$ in $\mathbf{Q}$.

Time Complexity Analysis: Three factors contribute to the time complexity of
AIV.GT. The first contributing factor is the learning of a PAG $\mathcal{P}$
from data and finding $\mathbf{S}$ from $\mathcal{P}$, which largely relies on
the rfci algorithm. In the worst situation, rfci has a complexity of
$\mathbf{O}(2^{r}*n)$, where $r$ is the maximum degree of a node in the
underlying causal MAG and $n$ is the sample size. In most cases, the average
degree of a causal Bayesian network is 2 to 5 [27], and most of the underlying
MAGs are sparse in real-world applications. Hence, the time complexity of rfci
is lower [28]. The complexity of Line 4 is $\mathbf{O}(1)$ since it reads from
$\mathcal{P}$. The second factor is estimating all possible causal effects,
i.e., Lines 8 to 11 in Algorithm 1. Noting that obtaining $\mathbf{Z}_{i}$
takes $\mathbf{O}(1)$ and calculating $\hat{\beta}_{i}$ needs
$\mathbf{O}(n*p^{2})$. Hence, the whole time complexity of this part is
$\mathbf{O}(|\mathbf{S}|*n*p^{2})$. The third factor is finding the pair of
IVs from $\mathbf{S}$ and time complexity relies on the size of $\mathbf{S}$
(pairwise search for a pair of VIs) and calculating covariance, which all
together takes $\mathbf{O}(|\mathbf{S}|^{2}*n*p^{2})$. Therefore, the overall
complexity of AIV.GT is
$\mathbf{O}(2^{r}*n+|\mathbf{S}|*n*p^{2}+|\mathbf{S}|^{2}*n*p^{2})=\mathbf{O}(2^{r}*n+|\mathbf{S}|^{2}*n*p^{2})$.
Therefore, the complexity of AIV.GT is largely attribute to the rfci algorithm
and searching for a pair satisfying the generalised tetrad criterion.

##  IV Experiments

We assess the performance of AIV.GT by comparing it to the state-of-the-art
causal effect estimators, firstly with a simulation study. Then, we conduct
experiments on two real-world datasets that have been used for a long time in
instrumental variable research [30, 31] to show that AIV.GT can be applied in
real-world applications.

Figure 2: The core of the true causal DAGs over measured and unmeasured used
to generate synthetic datasets. $\left|\mathbf{S}\right|$ denotes the number
of measured variables $Adj(W\cup Y)\setminus\\{W,Y\\}$. There are two valid
IVs $S_{1}$ and $S_{2}$ in all causal DAGs. In (a), $S_{1}$ and $S_{2}$ are
standard IVs. In (b), $S_{1}$ is a standard IV and $S_{2}$ is a CIV
conditioning on $X_{1}$. In (c), $S_{1}$ is a CIV conditioning on $\emptyset$
since $X_{1}$ is a collider and $S_{2}$ is a standard IV. In (d), $S_{1}$ and
$S_{2}$ are CIVs conditioning on $\emptyset$ and $X_{3}$, respectively. In
(e), $S_{1}$ is a condition IV conditioning on $\emptyset$ and $S_{2}$ is a
CIV conditioning on $\\{X_{3},X_{4}\\}$.

The estimators compared include (1) LSR, least squares regression $Y$ on
$\\{W,\mathbf{X}\\}$; (2) TSLS, two-stage least squares (TSLS [32]) of $Y$ on
$W$ using all variables $\mathbf{X}$ as standard IVs; (3) some invalid some
valid IV estimator (sisVIVE) [11]; (4) IV.tetrad method [8]. LSR is not an IV
based method. It is included since it is frequently used in Machine Learning
disregarding bias of latent variables in data, and it is used as a baseline.
All other three comparison estimators do not need a nominated IV. TSLS is a
standard IV estimator, and it is also used as a baseline. sisVIVE and
IV.tetrad are two most related methods and have been discussed in the
Introduction.

The Implementations of Estimators in Sections IV-A and IV-B2. The method OLS
is implemented by the function _cov_ in the $\mathbf{R}$ package _stats_. The
TSLS is programmed by the functions _cov_ in the $\mathbf{R}$ package _stats_
and _solve_ in the base. The implementation of LSR is same with TSLS, i.e.,
using the functions _cov_ and _solve_. The implementation of sisVIVE is based
on the function _sisVIVE_ in the $\mathbf{R}$ package _sisVIVE_. The
implementation of IV.tetrad is retrieved from the authors’
site111http://www.homepages.ucl.ac.uk/~ucgtrbd/code/iv_discovery. The
parameter of _num $\\_$ivs_ is set to 3 (2 for VitD). AIV.GT is implemented by
using the function _rfci_ in the $\mathbf{R}$ packages _pcalg_ , _cov_ in the
$\mathbf{R}$ package _stats_ , _solve_ in the base and the functions in
IV.tetrad.

The Implementations of Estimators in Section IV-B3. The implementations of the
compared estimators that require a known IV are introduced as follows. The
estimator TSLS is implemented by the function _ivreg_ in the $\mathbf{R}$
package _AER_ [33]. The implementation of TSLS.CIV is based on the functions
_glm_ and _ivglm_ from the $\mathbf{R}$ packages _stats_ and _ivtools_ [31].
FIVR is implemented by the function _instrumental_ $\\_$_forest_ in the
$\mathbf{R}$ package _grf_ [21]. All parameters in FIVR are default. AIViP is
implemented by the functions _rfci_ in $\mathbf{R}$ package _pcalg_ [26],
_glm_ in $\mathbf{R}$ package _stats_ and _ivglm_ in $\mathbf{R}$ package
_ivtools_ [31].

Evaluation Metrics & Parameter Setting. For the simulation study, we have the
ground truth of $\beta_{wy}$, so we report the estimation bias:
$\left|(\hat{\beta}_{wy}-\beta_{wy})/\beta_{wy}\right|*100$ (%). In the
experiments with real-world datasets, we empirically evaluate all estimators,
and then compare AIV.GT with four additional IV-based estimators that require
a nominated IV. The significant level $\alpha$ is set to 0.05 for the
functions of rfci and $Test.tetrad()$ in all experiments.

###  IV-A Simulation Study

The goal of this set of experiments is to test the effectiveness of AIV.GT
with and without colliders in the covariate set in comparison with various
causal effect estimators. We utilise five true DAGs over
$\mathbf{X}\cup\mathbf{U}\cup\\{W,Y\\}$ to generate five synthetic datasets
with latent variables. The five true DAGs are shown in Fig. 2. In addition to
the variables in the five true DAGs, 20 additional measured variables for each
dataset are generated as noise variables that are related to each other but
not to the variables in the five DAGs. The additional 20 variables are
generated by a multivariate normal distribution. Next, we will separately
introduce the details of each synthetic dataset generation based on the five
true DAGs in Fig. 2.

The synthetic dataset (a) is generated from the DAG (a) in Fig. 2, and the
specifications are as following: $U_{1}\sim Bernoulli(0.5)$ and
$S_{1},S_{2}\sim N(0,1)$, in which $N(,)$ denotes the normal distribution. The
treatment $W$ is generated from $n$ ($n$ denotes the sample size) Bernoulli
trials by using the assignment probability $P(W=1\mid
U_{1},S_{1},S_{2})=[1+exp\\{1-3*U_{1}-3*S_{1}-3*S_{2}\\}]$. The potential
outcome is generated from $Y_{W}=2+2*W+3*U_{1}+\epsilon_{w}$ where
$\epsilon_{w}\sim N(0,1)$.

The synthetic dataset (b) is generated from the DAG (b) in Fig. 2, and the
specifications are as following: $U_{1}\sim Bernoulli(0.5)$, $S_{1},S_{2}\sim
N(0,1)$ and $\epsilon_{X_{1}}\sim N(0,0.5)$. The measured variable
$X_{1}=0.8*S_{2}+\epsilon_{X_{1}}$.The treatment $W$ is generated by
$P(W=1\mid U_{1},S_{1},S_{2})=[1+exp\\{1-3*U_{1}-3*S_{1}-3*S_{2}\\}]$. The
potential outcome is generated from
$Y_{W}=2+2*W+3*U_{1}+3*X_{1}+\epsilon_{w}$.

The synthetic dataset (c) is generated from the DAG (c) in Fig. 2, and the
specifications are as following: $U_{1}\sim Bernoulli(0.5)$,
$S_{1},S_{2},U_{2},X_{2}\sim N(0,1)$ and $\epsilon_{X_{1}}\sim N(0,0.5)$. The
measured variable $X_{1}$ is generated by
$X_{1}=0.3+S_{1}+X_{2}+U_{2}+\epsilon_{X_{1}}$. The treatment $W$ is generated
by $P(W=1\mid U_{1},S_{1},S_{2})=[1+exp\\{1-3*U_{1}-3*S_{1}-3*S_{2}\\}]$. The
potential outcome is generated from
$Y_{W}=2+2*W+3*U_{1}+2*U_{2}+2*X_{2}+\epsilon_{w}$.

The synthetic dataset (d) is generated from the DAG (d) in Fig. 2, and the
specifications are as following: $U_{1}\sim Bernoulli(0.5)$,
$S_{1},S_{2},U_{2},X_{2}\sim N(0,1)$ and
$\epsilon_{X_{1}},\epsilon_{X_{3}}\sim N(0,0.5)$. The measured variables
$X_{1}$ and $X_{3}$ are generated by
$X_{1}=0.3+S_{1}+X_{2}+1.5*U_{2}+\epsilon_{X_{1}}$ and
$X_{3}=0.8*S_{2}+\epsilon_{X_{3}}$, respectively. The treatment $W$ is
generated by $P(W=1\mid
U_{1},S_{1},S_{2})=[1+exp\\{1-3*U_{1}-3*S_{1}-3*S_{2}\\}]$. The potential
outcome is generated by
$Y_{W}=2+2*W+3*U_{1}+2*U_{2}+2*X_{2}+2*X_{3}+\epsilon_{w}$.

The synthetic dataset (e) is generated from the DAG (e) in Fig. 2, and the
specifications are as following: $U_{1}\sim Bernoulli(0.5)$,
$S_{1},S_{2},U_{2},X_{2}\sim N(0,1)$ and
$\epsilon_{X_{1}},\epsilon_{X_{3}},\epsilon_{X_{4}}\sim N(0,0.5)$. The
measured variables $X_{1}$, $X_{3}$ and $X_{4}$ are generated by
$X_{1}=0.3+S_{1}+X_{2}+1.5*U_{2}+\epsilon_{X_{1}}$,
$X_{3}=0.8*S_{2}+\epsilon_{X_{3}}$ and $X_{4}=0.8*S_{2}+\epsilon_{X_{4}}$,
respectively. The treatment $W$ is generated by $P(W=1\mid
U_{1},S_{1},S_{2})=[1+exp\\{1-3*U_{1}-3*S_{1}-3*S_{2}\\}]$. The potential
outcome is generated by
$Y_{W}=2+2*W+3*U_{1}+2*U_{2}+2*X_{2}+2*X_{3}+2*X_{4}+\epsilon_{w}$.

The suitability of datasets with the method requirements is summarised as
Table I.

TABLE I: Satisfaction (tick)/violation (cross) of the assumptions of a method by a dataset. “$\checkmark$ ?” means the problem of collider bias suffered by IV.tetrad. | (a) | (b) | (c) | (d) | (e)
---|---|---|---|---|---
LSR | $\times$ | $\times$ | $\times$ | $\times$ | $\times$
TSLS | $\checkmark$ | $\times$ | $\times$ | $\times$ | $\times$
sisVIVE | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\times$ | $\times$
IV.tetrad | $\checkmark$ | $\checkmark$ |  $\checkmark$ ? |  $\checkmark$ ? |  $\checkmark$ ?
AIV.GT | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$

Figure 3: Estimation Bias (%) of the estimators on five synthetic datasets.
AIV.GT has the smallest bias on all datasets.

##### Results.

The estimation biases AIV.GT and the four compared estimators on the synthetic
datasets are visualised in Fig. 3. From Fig. 3, we have the following
observations: (1) The biases of LSR are large on all datasets. This is because
it does not consider any bias of latent variables in data. (2) TSLS has a low
bias on dataset (a) because $S_{1}$ and $S_{2}$ are standard IVs in this
dataset. It has large biases on all other datasets since there are no standard
IVs in the datasets. (3) sisVIVE works well on the first three datasets, i.e.
(a) $-$ (c), since its requirement (i.e. a half of covariates are IVs) is
satisfied. It works poorly on datasets (d) and (e) because of the assumption
is violated. (4) IV.tetrad has a low biases on the first two datasets, but
large biases on the last three datasets because the last three datasets
contain colliders in the covariate sets. All datasets satisfy the assumption
of IV.tetrad (i.e. a pair of CIVs), but it suffers the problem of collider
bias as identified in this paper. (5) AIV.GT obtains consistent results and
has the lowest bias on all datasets.

In sum, AIV.GT is able to obtain an unbiased $\hat{\beta}_{wy}$ from data with
latent variables when there exists a pair of AIVs. AIV.GT overcomes the
collider bias suffered by IV.tetrad.

###  IV-B Experiments on Two Real-world Datasets

It is challenging to evaluate the performance of causal effect estimators,
including AIV.GT on real-world datasets because the ground truth causal DAG
and $\beta_{wy}$ are not available. We select two real-world datasets with
empirical estimates available in the literature, including Vitamin D data
(VitD) [34, 31] and Schoolingreturns [30]. The two datasets have been
extensively studied and analysed before and each has a nominated AIV.
Therefore, it is credible to choose them as the benchmark datasets to evaluate
AIV.GT.

The two datasets have the nominated AIV wrt., $(W,Y)$, but the conditioning
sets are unknown. There is not an available algorithm in literature to
discover the conditioning set that instrumentalise the nominated AIV on both
datasets. Therefore, we divide the experiments on both datasets into two
parts: (1) experiments on AIV.GT in comparison with four estimators without
nominated AIVs; (2) experiments comparing AIV.GT with four additional IV
estimators that require a nominated AIV.

####  IV-B1 Details of the Two Real-world Datasets

##### Vitamin D (VitD).

This dataset was collected from a cohort study of vitamin D status on
mortality, i.e., the potential effect of VitD on death, reported in [34]. The
dataset contains 2571 individuals and 5 variables: age, filaggrin (a binary
variable indicating filaggrin mutations), vitd (a continuous variable measured
as serum 25-OH-D (nmol/L)), time (follow-up time), and death (binary outcome
indicating whether an individual died during follow-up) [31]. A measured value
of vitamin D less than 30 nmol/L implies vitamin D deficiency. We take the
estimated $\hat{\beta}_{wy}=2.01$ with 95% confidence interval $(0.96,4.26)$
from the literature [34] as the reference causal effect.

##### Schoolingreturns.

This dataset is from the national longitudinal survey of youth (NLSY) of US
young employees, aged range from 24 to 34 [30]. The dataset contains 3010
individuals and 19 variables. The treatment is the education of employees, and
the outcome is raw wages in 1976 (in cents per hour). The covariates include
experience (years of labour market experience), ethnicity (a factor indicating
ethnicity), resident information of an individual, age, nearcollege (whether
an individual grew up near a 4-year college), Education in 1966 (education66),
marital status, father’s educational attainment (feducation), mother’s
educational attainment (meducation), Ordered factor coding family education
class (fameducation), and so on. A goal of the studies on this dataset is to
investigate the causal effect of education on earnings. We take
$\hat{\beta}_{wy}=13.29$% with 95% confidence interval $(0.0484,0.2175)$ from
[35] as the reference causal effect.

####  IV-B2 Comparing AIV.GT with the Estimators without Requiring a Known IV

We conduct experiments on the two real-world datasets to assess AIV.GT against
the four estimators that do not require nominated IV as in Section IV-A with
the simulated data. All experimental results are visualised in Fig. 4 for VitD
and Schoolingreturns, respectively.

Figure 4: The experimental results of the five estimators without a given IV
for all estimators on two real-world datasets. The two dotted lines represent
empirically estimated causal effect with 95% confidence interval. Noting that
the estimated causal effect of LSR on VitD is close to zero and not visible in
the left panel.

##### Results on VitD.

From Fig. 4, we have the following observations: (1) the estimated result of
LSR is close to 0 and far away from the 95% confidence interval of the
empirical estimation; (2) the estimated results of TSLS, sisVIVE, IV.tetrad
and AIV.GT are close to the reference causal effect 2.01 and fall into the 95%
confidence interval.

##### Results on Schoolingreturns.

According to Fig. 4, we have the following findings: (1) the estimated result
of TSLS is at the bottom of the 95% confidence interval; (2) the estimated
results of LSR and sisVIVE fall outside of the empirical interval. It is very
likely that their assumptions have not been satisfied; (3) the estimated
results of IV.tetrad and AIV.GT are in the empirical interval. They are
consistency with the reference causal effect [35]. The consistency between the
results IV.tetrad and AIV.GT is likely due to the reason that they have found
proper conditional sets.

The experiments show that AIV-GT can obtain consistent estimations in both
real-world datasets.

####  IV-B3 Comparing AIV.GT with the Estimators with Known IVs

We add the four more comparison methods that require the given IVs, which are
(1) TSLS.IV, TSLS with a given IV; (2) TSLS.CIV [36], TSLS with a given CIV
$S$ by conditioning on $\mathbf{X}\setminus\\{S\\}$; (3) FIVR, causal random
forest for instrumental variable regression with a given CIV $S$ and
conditioning on $\mathbf{X}\setminus\\{S\\}$ [21]; (4) AIViP [24], Ancestral
IV estimator in PAG.

The two datasets have nominated IVs in the literature. The indicator of
filaggrin was used as an IV in VitD [34] and Card [30] used geographical
proximity to a college, i.e.,  _nearcollege_ as an IV in Schoolingreturns. All
results of the above four estimators and AIV.GT are visualised in Fig. 5.

Figure 5: The experimental results of five estimators with a given IV for all
comparison methods (AIV.GT does not use the given IV) on both real-world
datasets. The two dotted lines represent empirically interval with 95%
confidence interval.

##### Results on VitD.

AIV.GT discovers $\\{age,time\\}$ as a pair of AIVs. They are reasonable AIVs
since they affect $W$ (vitd, vitamin D status) but do not directly affect $Y$
(death). From Fig. 5, we see that the results of TSLS.CIV, AIViP and AIV.GT
are in the middle of 95% empirical interval of the reference causal effect.

##### Results on Schoolingreturns.

AIV.GT discovers $\\{feducation,fameducation\\}$ as a pair of AIVs. They are
valid IVs because father’s educational attainment and family education class
affect their child’s education, but do not directly affect the child’s income.
From Fig. 5, we observe that the results of TSLS.IV, AIViP and AIV.GT are in
the middle of the 95% empirical interval of the reference causal effect.

In a word, AIV.GT, which does not need given AIVs, performs better or
comparable with other methods which require a given IV. This shows the
potential of the AIV.GT in a broader range of real-world applications.

##  V Related Work

Latent variables are the major obstacle to estimating causal effect from
observational data [37, 6, 38]. When the treatment and outcome are confounded,
IV methods [5, 32, 39, 13, 6, 14, 40, 41] provide a solution.

Some methods have been developed for standard IV based causal effect
estimation when the IVs are given by domain experts [4, 39], such as the well-
known two-stage least squares IV estimator [32] (TSLS) which obtains causal
effect using the ratio of two regression coefficients. Recently, Athey et al.
[21] developed the generalised random forests to estimate conditional causal
effects by using non-parametric quantile regression and instrumental variable
regression (FIVR). The conditional causal effects can be aggregated into the
average causal effect and FIVR has been compared in our experiments. We refer
readers to [5, 23, 42] for a review of standard IV based methods.

When a CIV is given, a proper conditioning set needs to be identified for
unbiased causal effect estimation. Cheng et al.  [24] have proposed a data-
driven method AIViP for identifying such a conditioning set for casual effect
estimation with a given CIV. Our work is different from these works since we
focus on discovering AIVs and the corresponding conditioning set
simultaneously from data and the proposed method is more general than the
existing methods.

CIV approach is very similar to the approach of covariate adjustment since
both of them need to identify a proper conditioning set. However, there are
essential differences between CIV and covariate adjustment. Most methods for
covariate adjustment assume not a latent variable between $W$ on $Y$ [6, 43,
19]. There are four graphical criteria for identifying a proper conditioning
set from a causal graph: back-door criterion [6], adjustment criterion [44],
generalised back-door criterion [43] and generalised adjustment criterion
[19]. There are some data-driven methods based on the four graphical criteria
[45, 46, 47]. More detailed discussions for covariate selection can be found
in [48, 49, 50]. There are no works for finding a conditioning set when the
pair of $(W,Y)$ are confounded as discussed in this work.

##  VI Conclusion

Estimating causal effects in the presence of latent variables is a challenging
problem. IV is a well-known approach to address this challenge. However, most
existing IV methods require strong domain knowledge or assumptions to
determine an IV. This restricts the practical use of the IV approach. In this
paper, we present the theory and a practical algorithm (AIV.GT) for finding
valid AIVs and their corresponding conditioning sets from data, to enable
data-driven causal effect estimation from data with latent variables. The
experiments on synthetic datasets demonstrate that AIV.GT is able to address
the challenges of latent variables and outperform the state-of-the-art causal
effect estimators. The experimental results on two real-world datasets also
show that AIV.GT achieves consistent results with empirical estimates in the
literature, implying the practicability of AIV.GT in real-world applications.

## References

  * [1] O. Reiersøl, “Confluence analysis by means of instrumental sets of variables,” Ph.D. dissertation, Almqvist & Wiksell, 1945.
  * [2] G. W. Imbens and D. B. Rubin, _Causal Inference in Statistics, Social, and Biomedical Sciences_. Cambridge University Press, 2015.
  * [3] H. M. Blalock, _Causal Models in the Social Sciences_. Routledge, 2017.
  * [4] J. D. Angrist, G. W. Imbens, and D. B. Rubin, “Identification of causal effects using instrumental variables,” _Journal of the American Statistical Association_ , vol. 91, no. 434, pp. 444–455, 1996.
  * [5] E. P. Martens, W. R. Pestman _et al._ , “Instrumental variables: application and limitations,” _Epidemiology_ , vol. 17, no. 3, pp. 260–267, 2006.
  * [6] J. Pearl, _Causality_. Cambridge University Press, 2009.
  * [7] M. Kuroki and Z. Cai, “Instrumental variable tests for Directed Acyclic Graph Models,” in _International Conference on Artificial Intelligence and Statistics_ , pp. 190–197.
  * [8] R. Silva and S. Shimizu, “Learning instrumental variables with structural and non-gaussianity assumptions,” _Journal of Machine Learning Research_ , vol. 18, no. 120, pp. 1–49, 2017.
  * [9] J. Pearl, “On the testability of causal models with latent and instrumental variables,” in _Association for Uncertainty in Artificial Intelligence_ , 1995, pp. 435–443.
  * [10] T. Chu, R. Scheines, and P. Spirtes, “Semi-instrumental variables: a test for instrument admissibility,” in _Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence_ , 2001, pp. 83–90.
  * [11] H. Kang, A. Zhang _et al._ , “Instrumental variables estimation with some invalid instruments and its application to Mendelian randomization,” _Journal of the American Statistical Association_ , vol. 111, no. 513, pp. 132–144, 2016.
  * [12] J. S. Hartford, V. Veitch _et al._ , “Valid causal inference with (some) invalid instruments,” in _International Conference on Machine Learning_. PMLR, 2021, pp. 4096–4106.
  * [13] C. Brito and J. Pearl, “Generalized instrumental variables,” in _Proceedings of the Eighteenth Conference on Uncertainty in Artificial Intelligence_ , 2002, pp. 85–93.
  * [14] B. Van der Zander, M. Liśkiewicz, and J. Textor, “Efficiently finding conditional instruments for causal inference,” in _International Joint Conference on Artificial Intelligence_. AAAI Press, 2015, pp. 3243–3249.
  * [15] T. Richardson and P. Spirtes, “Ancestral graph markov models,” _The Annals of Statistics_ , vol. 30, no. 4, pp. 962–1030, 2002.
  * [16] J. Zhang, “Causal reasoning with ancestral graphs,” _Journal of Machine Learning Research_ , vol. 9, no. 7, pp. 1437–1474, 2008.
  * [17] P. Spirtes, C. N. Glymour _et al._ , _Causation, Prediction, and Search_. MIT Press, 2000.
  * [18] J. Zhang, “On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias,” _Artificial Intelligence_ , vol. 172, no. 16-17, pp. 1873–1896, 2008.
  * [19] E. Perković, J. Textor, and M. Kalisch, “Complete graphical characterization and construction of adjustment sets in markov equivalence classes of ancestral graphs,” _Journal of Machine Learning Research_ , vol. 18, pp. 1–62, 2018.
  * [20] A. Abadie, “Semiparametric instrumental variable estimation of treatment response models,” _Journal of Econometrics_ , vol. 113, no. 2, pp. 231–263, 2003.
  * [21] S. Athey, J. Tibshirani, and S. Wager, “Generalized random forests,” _The Annals of Statistics_ , vol. 47, no. 2, pp. 1148–1178, 2019.
  * [22] R. J. Bowden and D. A. Turkington, _Instrumental Variables_. Cambridge University Press, 1990, no. 8.
  * [23] M. A. Hernán and J. M. Robins, “Instruments for causal inference: an epidemiologist’s dream?” _Epidemiology_ , pp. 360–372, 2006.
  * [24] D. Cheng, J. Li _et al._ , “Ancestral instrument method for causal inference without complete knowledge,” in _International Joint Conference on Artificial Intelligence_ , 2022.
  * [25] C. F. Aliferis, A. Statnikov, I. Tsamardinos, S. Mani, and X. D. Koutsoukos, “Local causal and markov blanket induction for causal discovery and feature selection for classification part i: Algorithms and empirical evaluation,” _Journal of Machine Learning Research_ , vol. 11, no. Jan, pp. 171–234, 2010\.
  * [26] M. Kalisch, M. Mächler _et al._ , “Causal inference using graphical models with the R package pcalg,” _Journal of Statistical Software_ , vol. 47, no. 11, pp. 1–26, 2012.
  * [27] M. Scutari, “Learning Bayesian Networks with the bnlearn R package,” _Journal of Statistical Software_ , vol. 35, no. i03, pp. 1–22, 2010.
  * [28] D. Colombo, M. H. Maathuis _et al._ , “Learning high-dimensional directed acyclic graphs with latent and selection variables,” _The Annals of Statistics_ , vol. 40, no. 1, pp. 294–321, 2012.
  * [29] J. Wishart, “Sampling errors in the theory of two factors,” _British Journal of Psychology_ , vol. 19, no. 2, p. 180, 1928.
  * [30] D. Card, “Using geographic variation in college proximity to estimate the return to schooling,” 1993.
  * [31] A. Sjolander and T. Martinussen, “Instrumental variable estimation with the R package ivtools,” _Epidemiologic Methods_ , vol. 8, no. 1, 2019.
  * [32] J. D. Angrist and G. W. Imbens, “Two-stage least squares estimation of average causal effects in models with variable treatment intensity,” _Journal of the American statistical Association_ , vol. 90, no. 430, pp. 431–442, 1995\.
  * [33] W. H. Greene, _Econometric Analysis_. Pearson Education India, 2003.
  * [34] T. Martinussen _et al._ , “Instrumental variables estimation under a structural Cox model,” _Biostatistics_ , vol. 20, no. 1, pp. 65–79, 2019\.
  * [35] M. Verbeek, _A Guide to Modern Econometrics_. John Wiley & Sons, 2008.
  * [36] M. M. Glymour, E. J. Tchetgen Tchetgen, and J. M. Robins, “Credible mendelian randomization studies: approaches for evaluating the instrumental variable assumptions,” _American journal of epidemiology_ , vol. 175, no. 4, pp. 332–339, 2012.
  * [37] S. Greenland, “Quantifying biases in causal models: classical confounding vs collider-stratification bias,” _Epidemiology_ , pp. 300–306, 2003.
  * [38] W. Chen, R. Cai, K. Zhang, and Z. Hao, “Causal discovery in linear non-gaussian acyclic model with multiple latent confounders,” _IEEE Transactions on Neural Networks and Learning Systems_ , pp. 1–12, 2021.
  * [39] G. W. Imbens, “Instrumental Variables: An Econometrician’s Perspective,” _Statistical Science_ , vol. 29, no. 3, pp. 323–358, 2014\.
  * [40] N. Sokolovska and P.-H. Wuillemin, “The role of instrumental variables in causal inference based on independence of cause and mechanism,” _Entropy_ , vol. 23, no. 8, p. 928, 2021.
  * [41] F. Xie, Y. He _et al._ , “Testability of instrumental variables in linear non-gaussian acyclic causal models,” _Entropy_ , vol. 24, no. 4, p. 512, 2022\.
  * [42] M. Baiocchi, J. Cheng, and D. S. Small, “Instrumental variable methods for causal inference,” _Statistics in Medicine_ , vol. 33, no. 13, pp. 2297–2340, 2014.
  * [43] M. H. Maathuis, D. Colombo _et al._ , “A generalized back-door criterion,” _The Annals of Statistics_ , vol. 43, no. 3, pp. 1060–1088, 2015\.
  * [44] I. Shpitser, T. J. VanderWeele, and J. M. Robins, “On the validity of covariate adjustment for estimating causal effects,” in _International Conference on Uncertaintyin Artificial Intelligence_. AUAI Press, 2010, pp. 527–536.
  * [45] M. H. Maathuis, M. Kalisch, and P. Bühlmann, “Estimating high-dimensional intervention effects from observational data,” _The Annals of Statistics_ , vol. 37, no. 6A, pp. 3133–3164, 2009.
  * [46] D. Entner, P. Hoyer, and P. Spirtes, “Data-driven covariate selection for nonparametric estimation of causal effects,” in _Proceedings of Artificial Intelligence and Statistics_ , 2013, pp. 256–264.
  * [47] D. Cheng, J. Li _et al._ , “Toward unique and unbiased causal effect estimation from data with hidden variables,” _IEEE Transactions on Neural Networks and Learning Systems_ , 2022.
  * [48] B. C. Sauer, M. A. Brookhart, J. Roy, and T. VanderWeele, “A review of covariate selection for non-experimental comparative effectiveness research,” _Pharmacoepidemiology and Drug Safety_ , vol. 22, no. 11, pp. 1139–1145, 2013.
  * [49] J. Witte and V. Didelez, “Covariate selection strategies for causal inference: Classification and comparison,” _Biometrical Journal_ , vol. 61, no. 5, pp. 1270–1289, 2019.
  * [50] R. Guo, L. Cheng, J. Li, P. R. Hahn, and H. Liu, “A survey of learning causality with data: Problems and methods,” _ACM Computing Surveys (CSUR)_ , vol. 53, no. 4, pp. 1–37, 2020.