# The CMB Dipole: Eppur Si Muove

R. M. Sullivan∗ D. Scott  Physics and Astronomy, University of British
Columbia,
Vancouver, BC, Canada
∗E-mail<EMAIL_ADDRESS>
www.ubc.ca

###### Abstract

The largest temperature anisotropy in the cosmic microwave background (CMB) is
the dipole. The simplest interpretation of the dipole is that it is due to our
motion with respect to the rest frame of the CMB. As well as creating the
$\ell$=1 mode of the CMB sky, this motion affects all astrophysical
observations by modulating and aberrating sources across the sky. It can be
seen in galaxy clustering, and in principle its time derivative through a
dipole-shaped acceleration pattern in quasar positions. Additionally, the
dipole modulates the CMB temperature anisotropies with the same frequency
dependence as the thermal Sunyaev-Zeldovich (tSZ) effect and so these
modulated CMB anisotropies can be extracted from the tSZ maps produced by
Planck. Unfortunately this measurement cannot determine if the dipole is due
to our motion, but it does provide an independent measure of the dipole and a
validation of the $y$ maps. This measurement, and a description of the first-
order terms of the CMB dipole, are outlined here.

###### keywords:

Cosmic Microwave Background; Cosmic Microwave Background dipole; Special
relativity; Thermal Sunyaev-Zeldovich effect.

##  1 The CMB Sky from Planck

Planck111Planck (http://www.esa.int/Planck) is a project of the European Space
Agency (ESA) with instruments provided by two scientific consortia funded by
ESA member states and led by Principal Investigators from France and Italy,
telescope reflectors provided through a collaboration between ESA and a
scientific consortium led and funded by Denmark, and additional contributions
from NASA (USA). was a space-based telescope that measured the microwave sky
in nine wavebands, allowing it to capture not only the cosmic microwave
background (CMB) but also several Galactic and extragalactic foreground
components. This is most clearly seen in figure 51 from Ref. 1, which shows
the various wavebands of the Planck satellite and the frequency spectra of the
foreground signals across those bands. One signal of interest to this study is
the thermal Sunyaev-Zeldovich (tSZ) effect, which produces so-called $y$-type
distortion signals. This comes from CMB photons being inverse-Compton
scattered, mostly through hot galaxy clusters, which makes holes (or lowers
the flux) at low frequencies and up-scatters (or makes an excess flux) at high
frequencies. This signal allows us to construct a novel and independent
measure of the CMB dipole because temperature anisotropies stemming from the
CMB dipole contaminate the $y$ maps. It can also be used as a valuable test of
the quality of the $y$ maps. We will start in Sec. 2 by setting out relevant
notation for the unboosted CMB sky. Next, in Sec. 3 we will boost the CMB sky,
and explore the relevant terms that arise from that boost in the subsections.
Of particular relevance, in Sec. 3.3 we will discuss our measurement of the
dipole modulation terms that mix with the tSZ effect. We will finish in Sec. 4
with a short discussion and conclusion regarding our work.

##  2 The Unboosted CMB sky

To derive the connection between the $y$ map and the dipole we will begin by
defining some useful terms regarding the unboosted CMB sky:

$\displaystyle x$ $\displaystyle\equiv\frac{h\nu}{k_{\mathrm{B}}T};$ (1)
$\displaystyle I$
$\displaystyle\equiv\frac{2k^{3}_{\mathrm{B}}T^{3}}{h^{2}c^{2}}\frac{x^{3}}{e^{x}-1};$
(2) $\displaystyle f(x)$ $\displaystyle\equiv\frac{xe^{x}}{e^{x}-1};$ (3)
$\displaystyle Y(x)$ $\displaystyle\equiv x\frac{e^{x}+1}{e^{x}-1}-4.$ (4)

These are the dimensionless frequency, the standard Planck blackbody intensity
function, the frequency dependence of the CMB anisotropies and the relative
frequency dependence of the tSZ effect or $y$ type distortions, respectively.
Thus, to first order the anisotropies of intensity measured by Planck can be
written as

$\displaystyle\frac{\delta I(\hat{\mathbf{n}})}{If(x)}=\frac{\delta
T(\hat{\mathbf{n}})}{T_{\rm CMB}}+y(\hat{\mathbf{n}})Y(x),$ (5)

where $\hat{\mathbf{n}}$ is the line of sight direction on the sky and we have
only considered the temperature anisotropies and the $y$ signals here.

##  3 The Boosted CMB sky

If we apply a boost to Eq. 5, with a dimensionless velocity $\bm{\beta}$, we
find

$\displaystyle\frac{\delta I^{\prime}(\mathbf{\hat{n}^{\prime}})}{If(x)}=$
$\displaystyle\beta\mu+\frac{\delta T(\mathbf{\hat{n}^{\prime}})}{T_{\rm
CMB}}(1+3\beta\mu)$
$\displaystyle+Y(x)\left(y(\mathbf{\hat{n}^{\prime}})(1+3\beta\mu)+\beta\mu\frac{\delta
T(\mathbf{\hat{n}^{\prime}})}{T_{\rm CMB}}\right)$ $\displaystyle+\beta\mu
y(\mathbf{\hat{n}}^{\prime})\left(Y^{2}(x)-x\frac{dY(x)}{dx}\right)+\mathcal{O}(\beta^{2}),$
(6)

where $\mu=\cos(\theta)$, and $\theta$ is the angle between the boost
$\bm{\beta}$ and the line of sight $\mathbf{\hat{n}}^{\prime}$. The first line
has the same frequency dependence as thermal fluctuations and so appear in
typical CMB temperature anisotropy maps. Crucially for our analysis, the
middle line has the same frequency dependence as $y$-type distortions and thus
describes the signals in the $y$ map. The final line has more obscure
frequency dependence and is not discussed here. Additionally, the direction of
the incoming photons will change from $\mathbf{\hat{n}}$ to
$\mathbf{\hat{n}^{\prime}}$, where
$\mathbf{\hat{n}^{\prime}}=\mathbf{\hat{n}}-\nabla(\mathbf{\hat{n}}\cdot\mathbf{\beta})$;
this deflection of the photons by
$\nabla(\mathbf{\hat{n}}\cdot\mathbf{\beta})$ is due to aberration, an effect
that is not unique to the microwave sky and occurs for all astronomical
observations. We will now discuss each of these terms in turn.

###  3.1 The CMB dipole: $\beta\mu$

In the first line of Eq. 6, the term $\beta\mu$ describes the pure CMB dipole,
as discussed previously. This mainly (or perhaps entirely) comes from our
local motion with respect to the CMB rest frame and it has been previously
measured in Refs. 2, 3, and 4, and most recently in Refs. 5, 6, and 7. Taking
the large dipole as being solely caused by our motion, the velocity is
$v=(369.82\pm 0.11)$ km s-1 in the direction $(l,b)=(264$ .${}^{\circ}021\pm
0$ .${}^{\circ}011,48$ .${}^{\circ}253\pm 0$ .${}^{\circ}005)$ [5] and can be
easily seen in the CMB frequency maps, such as in Fig. 1.

Figure 1: Planck 100-GHz channel map from the NPIPE (PR4) data release [8],
showing the dominant $\ell=1$ mode or dipole across the sky. The temperature
difference across the sky here is 3.36 mK.

###  3.2 Aberration and Modulation of the CMB anisotropies:
$(1+3\beta\mu)\delta T(\mathbf{\hat{n}^{\prime}})/T_{\rm CMB}$

The second term in the first line of Eq. 6 is the dipole aberration and
modulation of the temperature anisotropies of the CMB. The modulation causes
the temperature anisotropies to be brighter in the forwards direction, and
dimmer in the reverse direction. The aberration causes the anisotropies to be
more condensed in the forwards direction, and more stretched out in the
reverse direction (effectively the same as $\ell=1$ lensing). These two
effects can be seen in Fig. 2. This effect was first measured in Ref. 9 to
about 4$\,\sigma$.

(a) The unboosted CMB sky.

(b) The modulated CMB sky.

(c) The aberrated CMB sky.

Figure 2: Here the CMB sky is shown unboosted in (a), with a modulation from a
boost of 90 % the speed of light in (b), and with aberration from a boost of
90 % of the speed of light in (c). In the case of modulation, the anisotropies
are more intense in the forward direction and less so in the reverse
direction, whereas the aberration condenses the anisotropies in the forwards
direction and causes them to be more spread out in the reverse direction.

###  3.3 Temperature modulation and the tSZ effect: $Y(x)\beta\mu\delta
T(\mathbf{\hat{n}^{\prime}})/T_{\rm CMB}$

The second line of Eq. 6 shows the dipole-generated signals in the $y$ maps
produced by Planck. The first half is the same modulation and aberration terms
as were seen in the temperature anisotropies; however, the final term is due
to the second-order expansion of the intensity about $T_{\rm CMB}$ and adds a
contribution to the $y$ maps from the temperature anisotropies.

We can look for this signal by cross-correlating a template map, derived from
the CMB temperature data, with a $y$ map. To this end, we use the so-called
2D-ILC CMB temperature map which was produced by the “Constrained ILC”
component-separation method designed by Ref. 10 to explicitly null out the
contribution from the $y$-type spectral distortions in the CMB map. We also
use the SMICA-NOSZ temperature map, similarly produced with the express intent
of removing the $y$-type spectral distortions, and which was generated for the
Planck 2018 data release [11]. Likewise, we use the corresponding 2D-ILC $y$
map, and the Planck MILCA $y$ map, which explicitly null out the contributions
from a (differential) blackbody spectral distribution in the $y$ map [12, 13].
If we multiply our CMB map with $\beta\mu$ and cross-correlate that with our
tSZ map, then we can directly probe the dipole modulation.

In Ref. 9 a quadratic estimator was used to determine the dipole aberration
and modulation, in essence using the auto-correlation of the CMB fluctuation
temperature maps. In this work we use the fact that we know the true CMB
fluctuations with excellent precision and therefore know the signal that
should be present in the $y$ map. Thus, we fully exploit the angular
dependence of the modulation signal and remove much of the cosmic variance
that would be present in the auto-correlation. In order to implement this idea
we define three templates, $B_{i}$ (with $i=1,2,3$) as

$\displaystyle B_{i}(\hat{\mathbf{n}})$
$\displaystyle=\beta\hat{\mathbf{n}}\cdot\hat{\mathbf{m}}_{i}\,\frac{\delta
T}{T_{0}}(\hat{\mathbf{n}}),$ (7)

where $\beta=v/c$ [5] is $1.23357\times 10^{-3}$ and
$\hat{\mathbf{m}}_{1},\hat{\mathbf{m}}_{2},\hat{\mathbf{m}}_{3}$ are the CMB
dipole direction, an orthogonal direction in the Galactic plane, and the third
remaining orthogonal direction (see Fig. 3). Due to the presence of the CMB
dipole, the signal $B_{1}$ should be present in the $y$ map and so we can
directly cross-correlate $B_{1}$ with our $y$ map to pull out the signal.
Likewise, the cross-correlation of $B_{2}$ and $B_{3}$ with our $y$ map should
give results consistent with noise.

Figure 3: Map of the tSZ effect from the MILCA component-separation method in
$y$-map units (top left) and the expected modulated CMB signal (top right)
generated using the SMICA-NOSZ CMB map in units of $T_{0}$. The bottom left
and right panels are the CMB anisotropies modulated in orthogonal directions
to the CMB dipole. Note that the map of the tSZ effect (top left) has a
different scale bar when compared to the other three (i.e., the modulation
signal is about 50 times weaker).

Our $y$ simulations are generated by first computing the power spectra of our
data $y$ maps; specifically we apply the MASTER method using the NaMASTER
routine [14] to account for the applied mask [15]. Then we generate $y$ maps
using this power-spectrum with the HEALPix [16] routine synfast. We finally
apply a Gaussian smoothing of 5′ to model the telescope beam. For each
analysis method we estimate the amplitude of the dipole ($\hat{\beta}_{i}$) in
each of the three orthogonal directions. We apply the same analysis procedure
on a suite of 1000 $y$ simulations, generated with and without the dipolar
modulation term.

We use two methods of cross-correlation: the first is performed directly in
map-space; and the second is performed in harmonic space. For both methods we
first apply our mask to the templates $B_{i}$ and the $y$ map.

In the map-space method we then locate all peaks (i.e., local maxima or
minima) of the template map $B_{i}$ and select a patch of radius $2$
.${}^{\circ}0$ around each peak. For every peak we obtain an estimate of
$\hat{\beta}_{i}$ through the simple operation

$\displaystyle\hat{\beta}_{i,p}$ $\displaystyle=\beta\frac{\sum_{k\in
D(p)}B_{i,k}y_{k}}{\sum_{k\in D(p)}B_{i,k}^{2}},$ (8)

where $D(p)$ is the collection of all _unmasked_ pixels in a $2$
.${}^{\circ}0$ radius centred on pixel $p$, and $p$ is the position of a peak.
Equation 8 is simply a cross-correlation in map space and by itself offers a
highly noisy (and largely unbiased) estimate. We then combine all individual
peak estimates with a set of weights ($w_{p}$) to give our full estimate:

$\displaystyle\hat{\beta}_{i}$
$\displaystyle=\frac{\sum_{p}w_{i,p}\hat{\beta}_{i,p}}{\sum_{p}w_{i,p}}.$ (9)

We choose $w_{p}$ to be proportional to the square of the dipole, and use
weights that are proportional to the square of the Laplacian at the peak [17];
this favours sharply defined peaks over shallow ones. Finally we account for
the scan strategy of the Planck mission by weighting by the 217-GHz hits
map[18], denoted $H^{217}_{p}$. The weights are then explicitly

$\displaystyle w_{i,p}$
$\displaystyle=|\hat{\mathbf{n}}\cdot\hat{\mathbf{m}}_{i}|^{2}_{p}\left(\left.\nabla^{2}(B_{i})\right|_{p}\right)^{2}H^{217}_{p}.$
(10)

We apply the method for each of our simulated $y$ maps, in exactly the same
way as for the data.

Under the assumption that the $y$ map contains the template ($B_{i}$), the $y$
multipoles are Gaussian random numbers with mean and variance given by

$\displaystyle s^{i}_{\ell m}$ $\displaystyle=\int
d\Omega\,\beta\,\hat{\mathbf{m}}_{i}\cdot\hat{\mathbf{n}}\,\frac{\delta
T}{T_{0}}M(\Omega)Y^{*}_{\ell m},$ (11) $\displaystyle\sigma^{2}_{\ell}$
$\displaystyle=C^{y}_{\ell}+N^{y}_{\ell},$ (12)

respectively, where $M(\Omega)$ is the mask over the sphere, $Y_{\ell m}$ are
the spherical harmonics, and the $\hat{\mathbf{m}}_{i}$ are as defined in Eq.
7. We can obtain an estimate of $\beta_{i}$ by taking the cross-correlation
with inverse-variance weighting. Our estimator is therefore

$\displaystyle\hat{\beta}_{i}$
$\displaystyle=\beta\sum_{i^{\prime}}\left[\sum_{\ell
m}^{\ell_{\max}}s^{i}_{\ell m}(s^{i^{\prime}}_{\ell
m})^{*}/\sigma^{2}_{\ell}\right]^{-1}\sum_{\ell
m}^{\ell_{\max}}s^{i^{\prime}}_{\ell m}(y_{\ell m})^{*}/\sigma^{2}_{\ell}.$
(13)

Figure 4: Histograms of $\hat{\beta}_{i}/\beta$ values (with 1, 2, and 3
corresponding to the CMB dipole direction, the Galactic plane, and a third
orthogonal direction) using the map-space analysis for MILCA $y$ maps, and for
CMB template maps SMICA-NOSZ. Blue histograms are simulations with the dipolar
modulation term, while orange histograms are simulations without modulation.
Black vertical lines denote the values of the real data, demonstrating that
they are much more consistent with the existence of the dipolar modulation
term than without it. Dashed lines show the 68 % regions for a Gaussian fit to
the histograms. To see the full results with all data analysis combinations
see Ref. 19. Figure 5: As in Fig. 4, except now for the harmonic-space
analysis.

First we compare the consistency of the data with our two sets of simulations
(with and without the dipole term). This comparisons shown in Figs. 4 and 5
have blue histograms being the simulations _with_ the dipole term and orange
histograms being _without_. The data (black line) for 2D-ILC and MILCA can
clearly be seen to be consistent with the simulations with the dipole term.
Further details and analysis may be found in Ref. 19.

##  4 Conclusion: Distinguishing Intrinsic and Extrinsic CMB Dipoles

The frequency-dependent part of the dipolar-modulation signal is agnostic to
the source of the large CMB dipole. Therefore, its measurement is an
independent determination of the CMB dipole. While it may be tempting to use
this signal to distinguish an intrinsic dipole, it has been shown that an
intrinsic dipole and a dipole induced by a velocity boost would in fact have
the same dipolar-modulation signature on the sky [20, 21]

Due to the existence of the CMB dipole, a tSZ map necessarily contains a
contaminating signal that is simply the dipole modulation of the CMB
anisotropies. This occurs because CMB experiments do not directly measure
temperature anisotropies, but instead measure intensity variations that are
conventionally converted to temperature variations. This contamination adds
power to the tSZ map in a $Y_{20}$ pattern, with its axis parallel to the
dipole direction. We have measured this effect and determined a statistically
independent value of the CMB dipole, which is consistent with direct
measurements of the dipole. Using a conservative multipole cut on the $y$ map,
the significance of the detection of the dipole modulation signal is around 5
or $6\,\sigma$, depending on the precise choice of data set and analysis
method.

The question as to whether an intrinsic dipole could ever be observationally
distinguished from an extrinsic dipole (i.e. a Doppler boost) remains an open
question. The terms discussed in Eq. 6 are based on the assumption of a CMB
blackbody spectrum and cannot be used to distinguish the two, as they would
naturally arise whether the CMB dipole is caused by a boost, or if there is
for some other reason a dipole on the sky with the same magnitude and
direction.

## Acknowledgements

We would like to acknowledge the support of the Natural Sciences and
Engineering Research Council of Canada. Some of the results in this paper have
been derived using the HEALPix package. Results are based on observations
obtained with Planck (http://www.esa.int/Planck), an ESA science mission with
instruments and contributions directly funded by ESA Member States, NASA, and
Canada. We would also like to thank Dagoberto Contreras for collaboration on
topics within this paper.

## References

  * [1] Planck Collaboration X, Planck 2015 results. X. Diffuse component separation: Foreground maps, A&A 594, p. A10 (2016).
  * [2] A. Kogut, C. Lineweaver, G. F. Smoot, C. L. Bennett, A. Banday, N. W. Boggess, E. S. Cheng, G. de Amici, D. J. Fixsen, G. Hinshaw, P. D. Jackson, M. Janssen, P. Keegstra, K. Loewenstein, P. Lubin, J. C. Mather, L. Tenorio, R. Weiss, D. T. Wilkinson and E. L. Wright, Dipole Anisotropy in the COBE Differential Microwave Radiometers First-Year Sky Maps, ApJ 419, p. 1 (December 1993).
  * [3] D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer and E. L. Wright, The Cosmic Microwave Background Spectrum from the Full COBE FIRAS Data Set, ApJ 473, p. 576 (December 1996).
  * [4] G. Hinshaw, J. L. Weiland, R. S. Hill, N. Odegard, D. Larson, C. L. Bennett, J. Dunkley, B. Gold, M. R. Greason, N. Jarosik, E. Komatsu, M. R. Nolta, L. Page, D. N. Spergel, E. Wollack, M. Halpern, A. Kogut, M. Limon, S. S. Meyer, G. S. Tucker and E. L. Wright, Five-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Data Processing, Sky Maps, and Basic Results, ApJS 180, 225 (February 2009).
  * [5] Planck Collaboration I, Planck 2018 results. I. Overview, and the cosmological legacy of Planck, A&A 641, p. A1 (2020).
  * [6] Planck Collaboration II, Planck 2018 results. II. Low Frequency Instrument data processing, A&A, in press (2019).
  * [7] Planck Collaboration III, Planck 2018 results. III. High Frequency Instrument data processing, A&A, in press (2019).
  * [8] Planck Collaboration Int. LVII, Planck intermediate results. LVII. NPIPE: Joint Planck LFI and HFI data processing, A&A 643, p. 42 (2020).
  * [9] Planck Collaboration XXVII, Planck 2013 results. XXVII. Doppler boosting of the CMB: Eppur si muove, A&A 571, p. A27 (2014).
  * [10] M. Remazeilles, J. Delabrouille and J.-F. Cardoso, CMB and SZ effect separation with constrained Internal Linear Combinations, MNRAS 410, 2481 (February 2011).
  * [11] Planck Collaboration IV, Planck 2018 results. IV. Diffuse component separation, A&A, in press (2019).
  * [12] G. Hurier, J. F. Macías-Pérez and S. Hildebrandt, MILCA, a modified internal linear combination algorithm to extract astrophysical emissions from multifrequency sky maps, A&A 558, p. A118 (October 2013).
  * [13] Planck Collaboration XXII, Planck 2015 results. XXII. A map of the thermal Sunyaev-Zeldovich effect, A&A 594, p. A22 (2016).
  * [14] D. Alonso, J. Sanchez, A. Slosar and L. D. E. S. Collaboration, A unified pseudo-C$\ell$ framework, Monthly Notices of the Royal Astronomical Society 484, 4127 (01 2019).
  * [15] E. Hivon, K. M. Górski, C. B. Netterfield, B. P. Crill, S. Prunet and F. Hansen, MASTER of the Cosmic Microwave Background Anisotropy Power Spectrum: A Fast Method for Statistical Analysis of Large and Complex Cosmic Microwave Background Data Sets, ApJ 567, 2 (March 2002).
  * [16] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke and M. Bartelmann, HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere, ApJ 622, 759 (April 2005).
  * [17] V. Desjacques, Baryon acoustic signature in the clustering of density maxima, Phys. Rev. D 78, p. 103503 (November 2008).
  * [18] Planck Collaboration VIII, Planck 2015 results. VIII. High Frequency Instrument data processing: Calibration and maps, A&A 594, p. A8 (2016).
  * [19] Planck Collaboration Int. LVI, Planck intermediate results. LVI. Detection of the CMB dipole through modulation of the thermal Sunyaev-Zeldovich effect: Eppur si muove II, A&A 644, p. 100 (2020).
  * [20] A. Challinor and F. van Leeuwen, Peculiar velocity effects in high-resolution microwave background experiments, Phys. Rev. D 65, p. 103001 (May 2002).
  * [21] A. Notari and M. Quartin, CMB all-scale blackbody distortions induced by linearizing temperature, ArXiv e-prints (October 2015).


# Identifying Counterfactual Queries with the R Package cfid

Santtu Tikka
University of Jyväskylä<EMAIL_ADDRESS>
Santtu Tikka Identifying Counterfactual Queries with the R Package cfid
Identifying Counterfactual Queries In the framework of structural causal
models, counterfactual queries describe events that concern multiple
alternative states of the system under study. Counterfactual queries often
take the form of “what if” type questions such as “would an applicant have
been hired if they had over 10 years of experience, when in reality they only
had 5 years of experience?” Such questions and counterfactual inference in
general is crucial, for example when addressing the problem of fairness in
decision-making. Because counterfactual events contain contradictory states of
the world, it is impossible to conduct a randomized experiment to address them
without making several restrictive assumptions. However, it is sometimes
possible to identify such queries from observational and experimental data by
representing the system under study as a causal model, and the available data
as symbolic probability distributions. Shpitser and Pearl (2007) constructed
two algorithms, called ID* and IDC*, for identifying counterfactual queries
and conditional counterfactual queries, respectively. These two algorithms are
analogous to the ID and IDC algorithms by Shpitser and Pearl (2006b, a) for
identification of interventional distributions, which were implemented in R by
Tikka and Karvanen (2017) in the causaleffect package. We present the R
package cfid that implements the ID* and IDC* algorithms. Identification of
counterfactual queries and the features of cfid are demonstrated via examples.
causality, causal model, counterfactual, do-calculus, graph, identifiability
causality, causal model, counterfactual, do-calculus, graph, identifiability
Santtu Tikka
Department of Mathematics Statistics
Faculty of Mathematics and Science
University of Jyväskylä
P.O.Box 35, FI-40014, Finland
E-mail:
URL: http://users.jyu.fi/~santikka/

##  1 Introduction

Pearl’s ladder of causation (or causal hierarchy) consists of three levels:
association, intervention, and counterfactual (Pearl, 2009). These levels
describe a hierarchy of problems in increasing conceptual and formal
difficulty. On the first and lowest level, inference on associations is based
entirely on observed data in the form of questions such as “what is the
probability that an event occurs?” or “what is the correlation between two
variables”. On the second level, the inference problems are related to
manipulations of the system under study such as “what is the probability of an
event if we change the value of one variable in the system”. Questions on the
intervention level cannot be answered using tools of the association level,
because simply observing a change in a system is not the same as intervening
on the system. Randomized controlled trials are the gold standard for studying
the effects of interventions, because they enable the researcher to account
for confounding factors between the treatment and the outcome and to carry out
the intervention in practice. However, there are often practical limitations
that make it difficult, expensive, or impossible to conduct a randomized
experiment. The third and highest level is the counterfactual level.
Typically, counterfactual statements compare the real world, where an action
was taken or some event was observed, to an alternative hypothetical scenario,
where a possibly different action was taken, or a different event was
observed. Counterfactuals are often challenging to understand even
conceptually due this notion of contradictory events in alternative worlds,
and such alternatives need not be limited to only two. In general, questions
on the counterfactual level cannot be answered by relying solely on the
previous levels: no intervention or association is able to capture the notion
of alternative hypothetical worlds.

While counterfactual statements can be challenging, they are a core part of
our everyday thinking and discourse. Importantly, counterfactuals often
consider retrospective questions about the state of the world, such as “would
an applicant have been hired if they had more work experience”. This kind of
retrospection is crucial when fair treatment of individuals is considered in
hiring, healthcare, receiving loans or insurance, etc., with regards to
protected attributes, especially when the goal is automated decision-making.
Statistical approaches to fairness are insufficient in most contexts, such as
in scenarios analogous to the well-known Simpson’s paradox, but routinely
resolved using the framework of causal inference. In some cases, even
interventional notions of fairness may be insufficient, necessitating
counterfactual fairness (Kusner _et al._ , 2017; Zhang and Bareinboim, 2018).

The structural causal model (SCM) framework of Pearl provides a formal
approach to causal inference of interventional and counterfactual causal
queries (Pearl, 2009). An SCM represents the system of interest in two ways,
First, the causal relationships are depicted by a directed acyclic graph (DAG)
whose vertices correspond to variables under study and whose edges depict the
direct functional causal relationships between the variables. Typically, only
some of these variables are observed and the remaining variables are
considered latent, corresponding either to confounders between multiple
variables or individual random errors of single variables. Second, the
uncertainty related to the variables in the system is captured by assuming a
joint probability distribution over its latent variables. The functional
relationships of the model induce a joint probability distribution over the
observed variables. The SCM framework also incorporates the notion of external
interventions symbolically via the do-operator, and a graphical representation
of counterfactual scenarios via parallel worlds graphs (Avin _et al._ , 2005;
Shpitser and Pearl, 2007, 2008).

One of the fundamental problems of causal inference is the so-called
identifiability problem, especially the identifiability of interventional
distributions. Using the SCM framework and do-calculus, it is sometimes
possible to uniquely represent an interventional distribution using only the
observed joint probability distribution of the model before the intervention
took place. Such interventional distributions are called _identifiable_. More
generally, we say that a causal query is identifiable, if it can be uniquely
represented using the available data. In most identifiability problems, the
available data consists of causal quantities on levels below the query in the
ladder of causation, but the levels also sometimes overlap, (e.g., Bareinboim
and Pearl, 2012; Tikka and Karvanen, 2019; Lee _et al._ , 2019). The
identifiability problem of interventional distributions, and many other
interventional identifiability problems have been solved by providing a sound
and complete identification algorithm (e.g., Shpitser and Pearl, 2006b; Huang
and Valtorta, 2006; Lee _et al._ , 2019; Kivva _et al._ , 2022).

Software for causal inference is becoming increasingly prominent. For R (R
Core Team, 2022), a comprehensive overview of the state-of-the-art is provided
by the recently launched task view on causal inference
(https://cran.r-project.org/web/views/CausalInference.html) on the
Comprehensive R Archive Network (CRAN). Out of the packages listed in this
task view, the Counterfactual (Chen _et al._ , 2020) and WhatIf (Stoll _et
al._ , 2020) packages are directly linked to counterfactual inference, but the
focus of these packages is estimation and they do not consider the
identifiability of counterfactual queries. The R6causal (Karvanen, 2022)
package can be used to simulate data from counterfactual scenarios in a causal
model. R packages most closely related to causal identifiability problems are
the causaleffect (Tikka and Karvanen, 2017), dosearch (Tikka _et al._ , 2021),
and dagitty (Textor _et al._ , 2017). In Python (Van Rossum and Drake, 2009),
the Ananke module (Nabi _et al._ , 2020; Lee and Shpitser, 2020; Bhattacharya
_et al._ , 2020) and the DoWhy library (Sharma _et al._ , 2019) provide
comprehensive tools for causal inference. However, all the aforementioned
packages perform identification at the intervention level, not counterfactual.

We present the first implementation of the counterfactual identifiability
algorithms of Shpitser and Pearl (2007) (see also Shpitser and Pearl, 2008) as
the R package cfid (counterfactual identification). The cfid package also
provides a user-friendly interface for defining causal diagrams and the
package is compatible for other major R packages for causal identifiability
problems such as causaleffect, dosearch and dagitty by supporting graph
formats used by these packages as inputs.

The paper is organized as follows. Section 2 introduces the notation, core
concepts and definitions, and provides an example on manual identification of
a counterfactual query without relying on the identifiability algorithms.
Section 3 presents the algorithms implemented in cfid and demonstrates their
functionality via examples by tracing their operation line by line. Section 4
demonstrates the usage of the cfid package in practice. Section 5 concludes
the paper with a summary.

##  2 Notation and definitions

We follow the notation used by Shpitser and Pearl (2008) and we assume the
reader to be familiar with standard graph theoretic concepts such as ancestral
relations between vertices and d-separation. We use capital letters to denote
random variables and lower-case letters to denote their value assignments.
Bold letters are used to denote sets of random variables and counterfactual
variables. We associate the vertices of graphs with their respective random
variables and value assignments in the underlying causal models. In figures,
observed variables of graphs are denoted by circles, variables fixed by
interventions are denoted by squares, and latent unobserved variables are
denoted by dashed circles when explicitly included and by bidirected edges
when the corresponding latent variable has two observed children. Latent
variables with only one child, which are called _error terms_ , are not shown
for clarity.

A _structural causal model_ is a tuple
$M=(\mathbf{U},\mathbf{V},\mathbf{F},P(\mathbf{u}))$ where $\mathbf{U}$ is a
set of unobserved random variables, $\mathbf{V}$ is a set of $n$ observed
random variables, $\mathbf{F}$ is a set of $n$ functions such that each
function $f_{i}$ is a mapping from
$\mathbf{U}\cup\mathbf{V}\setminus\\{V_{i}\\}$ to $V_{i}$ and such that it is
possible to represent the set $\mathbf{V}$ as function of $\mathbf{U}$.
$P(\mathbf{u})$ is a joint probability distribution over $\mathbf{U}$. The
causal model also defines its causal diagram $G$. Each $V_{i}\in\mathbf{V}$
corresponds to a vertex in $G$, and there is a directed edge from each
$V_{j}\in\mathbf{U}\cup\mathbf{V}\setminus\\{V_{i}\\}$ to $V_{i}$. We restrict
our attention to _recursive_ causal models in this paper, meaning models that
induce an acyclic causal diagram.

A _counterfactual variable_ $Y_{\mathbf{x}}$ denotes the variable $Y$ in the
submodel $M_{\mathbf{x}}$ obtained from $M$ by forcing the random variables
$\mathbf{X}$ to take the values $\mathbf{x}$ (often denoted by the do-operator
as $\textrm{do}(\mathbf{X}=\mathbf{x})$ or simply $\textrm{do}(\mathbf{x})$).
The distribution of $\mathbf{Y}_{\mathbf{x}}$ in the submodel $M_{\mathbf{x}}$
is called the _interventional distribution_ of $\mathbf{Y}$ and it is denoted
by $P_{\mathbf{x}}(\mathbf{y})$. However, if we wish to consider multiple
counterfactual variables that originate from different interventions, we must
extend our notation to counterfactual conjunctions. _Counterfactual
conjunctions_ are constructed from value assignments of counterfactual
variables, and individual assignments are separated by the $\wedge$ symbol.
For example, $y_{x}\wedge z_{x}\wedge x^{\prime}$ denotes the event that
$Y_{x}=y$, $Z_{x}=z$ and $X=x^{\prime}$. The probability $P(y_{x}\wedge
z_{x}\wedge x^{\prime})$ is the probability of the counterfactual event. Note
that primes do not differentiate variables, instead they are used to
differentiate between values i.e., $x$ is a different value from $x^{\prime}$
and they are both different from $x^{\prime\prime}$ but all three are value
assignments of the random variable $X$. If the subscript of each variable in
the conjunction is the same, the counterfactual probability simply reduces to
an interventional distribution.

Each counterfactual conjunction is associated with multiple _parallel worlds_
, each induced by a unique combination of subscripts that appears in the
conjunction. A _parallel worlds graph_ of the conjunction is obtained by
combining the graphs of the submodels induced by interventions such that the
latent variables are shared. The simplest version of a parallel worlds graph
is a twin network graph, contrasting two alternative worlds (Balke and Pearl,
1994a, b; Avin _et al._ , 2005). As a more complicated example, consider the
counterfactual conjunction $\gamma=y_{x}\wedge x^{\prime}\wedge z_{d}\wedge
d$. In simpler terms, this conjunction states that $Y$ takes the value $y$
under the intervention $\textrm{do}(X=x)$, $Z$ takes the value $z$ under the
intervention $\textrm{do}(D=d)$, and $X$ and $D$ take the values $x^{\prime}$
and $d$, respectively, when no intervention took place. Importantly, this
conjunction induces three distinct parallel worlds: the non-interventional (or
observed) world, a world where $X$ was intervened on, and a world where $D$
was intervened on. For instance, if the graph in Figure 1(a) depicts the
original causal model over the variables $Y,X,Z,W$ and $D$, then Figure 1(b)
shows the corresponding parallel worlds graph for $\gamma$, where each
distinct world is represented by its own set of copies of the original
variables. In Figure 1(b), $U$ corresponds to the bidirected edge between $X$
and $Y$ in Figure 1(a), and the other $U$-variables are the individual error
terms of each observed variable, that are not drawn when they have only one
child in Figure 1(a).

Note that instead of random variables, some nodes in the parallel worlds graph
now depict fixed values as assigned by the interventions in the conjunction.
This is a crucial aspect when d-separation statements are considered between
counterfactual variables in the parallel worlds graph, as a backdoor path
through a fixed value is not open. Furthermore, not every variable is
necessarily unique in a parallel worlds graph, making it possible to obtain
misleading results if d-separation is used to infer conditional independence
relations between counterfactual variables. For instance, if we consider the
counterfactual variables $Y_{x}$, $D_{x}$ and $Z$ in a causal model whose
diagram is the graph shown in Figure 1(a), then $Y_{x}$ is independent of
$D_{x}$ given $Z$, even though $Y_{x}$ is not d-separated from $D_{x}$ in the
corresponding parallel worlds graph of Figure 1(b). This conditional
independence holds because $Z$ and $Z_{x}$ are in fact the same counterfactual
variable. To overcome this problem, the parallel worlds graph must be further
refined into the _counterfactual graph_ where every variable is unique, which
we will discuss in the following sections in more detail. For causal diagrams
and counterfactual graphs, $V(G)$ denotes the set of observable random
variables not fixed by interventions and $v(G)$ denotes the corresponding set
of value assignments.

0$\scriptstyle X$0$\scriptstyle W$0$\scriptstyle D$0$\scriptstyle
Z$0$\scriptstyle Y$0 (a) A causal diagram.

0$\scriptstyle X$0$\scriptstyle W$0$\scriptstyle D$0$\scriptstyle
Z$0$\scriptstyle Y$0$x$0$\scriptstyle W_{x}$0$\scriptstyle
D_{x}$0$\scriptstyle Z_{x}$0$\scriptstyle Y_{x}$0$\scriptstyle
X_{d}$0$\scriptstyle W_{d}$0$d$0$\scriptstyle Z_{d}$0$\scriptstyle
Y_{d}$0$\scriptstyle U_{D}$0$\scriptstyle U$0$\scriptstyle U$0$\scriptstyle
U_{Z}$0$\scriptstyle U_{W}$ (b) A parallel worlds graph of (a) for
$y_{x}\wedge x^{\prime}\wedge z_{d}\wedge d$. Colors are used here to
distinguish the observed and fixed nodes that belong to different parallel
worlds: black for the non-interventional world, blue for the world induced by
$\textrm{do}(X=x)$, and red for the world induced by $\textrm{do}(D=d)$. Note
that node $U$ is drawn twice for clarity due to its many endpoints.

Figure 1: An example causal diagram and a corresponding parallel worlds graph.

The following operations are defined for counterfactual conjunctions and sets
of counterfactual variables: $\mathrm{sub}(\cdot)$ returns the set of
subscripts, $\mathrm{var}(\cdot)$ returns the set of (non-counterfactual)
variables, and $\mathrm{ev}(\cdot)$ returns the set of values (either fixed by
intervention or observed). For example, consider again the conjunction
$\gamma=y_{x}\wedge x^{\prime}\wedge z_{d}\wedge d$. Now,
$\mathrm{sub}(\gamma)=\\{x,d\\}$, $\mathrm{var}(\gamma)=\\{Y,X,Z,D\\}$ and
$\mathrm{ev}(\gamma)=\\{y,x,x^{\prime},z,d\\}$. Finally, $\mathrm{val}(\cdot)$
is the value assigned to a given counterfactual variable, e.g.,
$\mathrm{val}(y_{x})=y$. The notation $y_{\mathbf{x}..}$ denotes a
counterfactual variable derived from $Y$ with the value assignment $y$ in a
submodel $M_{\mathbf{x}\cup\mathbf{z}}$ where
$\mathbf{Z}\subseteq\mathbf{V}\setminus\mathbf{X}$ is arbitrary.

The symbol $P_{*}$ is used to denote the set of all interventional
distributions of a causal model $M$ over a set of observed variables
$\mathbf{V}$, i.e.,

$P_{*}=\\{P_{\mathbf{x}}\mid\mathbf{x}\text{ is any value assignment of
}\mathbf{X}\subseteq\mathbf{V}\\}$

In the following sections, we consider identifiability of counterfactual
queries in terms of $P_{*}$. In essence, this means that a counterfactual
probability distribution $P(\gamma)$ is identifiable if it can be expressed
using purely interventional and observational probabilities of the given
causal model.

###  2.1 Example on identifiability of a counterfactual query

We consider the identifiability of the conditional counterfactual query
$P(y_{x}|z_{x}\wedge x^{\prime})$ from $P_{*}$ in the graph depicted in Figure
2. This graph could for instance depict the effect of an applicant’s education
($X$) on work experience ($Z$) and a potential hiring decision ($Y$) by a
company. Our counterfactual query could then consider the statement “what is
the probability to be hired if the applicant’s education level was changed to
$x$, given that their work experience under the same intervention was $z$ and
when in reality their education level was $x^{\prime}$”. In this example, we
will not rely on any identifiability algorithms. Instead, we can derive a
formula for the counterfactual query as follows:

0$\scriptstyle X$0$\scriptstyle Z$0$\scriptstyle Y$ Figure 2: A graph for the
example on identifiability of a conditional counterfactual query
$P(y_{x}|z_{x}\wedge x^{\prime})$.

$\displaystyle P(y_{x}|z_{x}\wedge x^{\prime})$
$\displaystyle=\frac{P(y_{x}\wedge z_{x}\wedge
x^{\prime})}{\sum_{y}P(y_{x}\wedge z_{x}\wedge x^{\prime})}$
$\displaystyle=\frac{P(y_{xz}\wedge z_{x}\wedge
x^{\prime})}{\sum_{y}P(y_{xz}\wedge z_{x}\wedge x^{\prime})}$
$\displaystyle(\text{composition})$ $\displaystyle=\frac{P(y_{xz}|z_{x}\wedge
x^{\prime})P(z_{x}\wedge x^{\prime})}{\sum_{y}P(y_{xz}|z_{x}\wedge
x^{\prime})P(z_{x}\wedge x^{\prime})}$
$\displaystyle=\frac{P(y_{xz})P(z_{x}\wedge x^{\prime})}{P(z_{x}\wedge
x^{\prime})\sum_{y}P(y_{xz})}$ $\displaystyle(\text{independence
restrictions})$ $\displaystyle=P(y_{xz})$ $\displaystyle=P_{xz}(y)$

Thus we find that the answer to our initial question is simply the probability
of hiring if the applicant’s education level and work experience were changed
to $x$ and $z$, respectively. In the above derivation, we used the notions of
composition and independence restrictions (Holland, 1986; Pearl, 1995;
Halpern, 1998; Pearl, 2009). Composition is one of the axioms of
counterfactuals stating that if a variable is forced to a value that it would
have taken without the intervention, then the intervention will not affect
other variables in the system. In this case, intervention setting $Z_{x}$ to
$z$ has no effect on $Y_{x}$ because we have observed $Z_{x}=z$, thus we can
add $Z$ to the intervention set of $Y_{x}$. Independence restrictions state if
the observed parents of a variable are intervened on, then the counterfactual
is independent of any other observed variable when their parents are also held
fixed, if there are no paths between the variables via latent variables. In
this case $Y_{x,z}$ is independent of $Z_{x}$ and $X$ because there is no path
via latent variables connecting $Y$ to $Z$ or $X$ in $G$.

In this example, the interventional distribution $P_{x,z}(y)$ can be further
identified from the observed joint distribution $P(x,z,y)$ as $P(y|x,z)$ via
the second rule of do-calculus by noting that $Y$ is d-separated from $X$ and
$Z$ in the graph when the outgoing edges of $X$ and $Z$ are removed. Thus the
answer to our initial question can be further refined into the probability of
hiring among applicants with education level $x$ and work experience $z$. The
cfid package provides this kind of identification pipeline from the
counterfactual level down to the lowest possible level in the causal
hierarchy.

##  3 Algorithms for identifying counterfactual queries

Manual identification of counterfactuals is challenging and more nuanced than
identification of interventional distributions due to fixed values and non-
uniqueness of counterfactual variables in the parallel worlds graph.
Therefore, identification of a counterfactual query can be achieved in several
ways. First, we may find that the query is identifiable and thus we can
express it in terms of purely interventional distributions. In contrast, we
may find that the query is not identifiable, meaning that is not possible to
represent it in terms of purely interventional distributions. Alternatively,
we may find that the query is _inconsistent_ meaning that the same
counterfactual variable has been assigned at least two different values in the
conjunction, and thus the query is identified as a zero-probability event. For
example, suppose we are tasked with identifying $P(y_{x},y^{\prime}_{z})$ but
we find that $Y_{x}$ and $Y_{z}$ are actually the same variable, and thus
cannot attain two different values $y$ and $y^{\prime}$ simultaneously. For
conditional counterfactual queries, there is also a fourth option where the
query is undefined if the conditioning conjunction is inconsistent.

Algorithmic identification of interventional distributions takes advantage of
the so-called _C-component factorization_ (Tian and Pearl, 2002; Shpitser and
Pearl, 2006b) which also plays a key role in the identification of
counterfactual queries. The _maximal C-components_ of a causal diagram are
obtained by partitioning the vertices $\mathbf{V}$ related to observed
variables of the graph such that two vertices $A,B\in B$ in the same partition
are connected by a path with edges into $A$ and $B$ where every node on the
path in $\mathbf{V}$ except $A$ and $B$ is a collider, and $A$ and $B$ are not
connected to any other partitions via such paths. Maximal C-components are
defined analogously for parallel worlds graphs and counterfactual graphs with
the restriction that we do not consider vertices that correspond to fixed
values to belong to any C-component. The set of maximal C-components of a DAG
$G$ is denoted by $C(G)$. As an example, the maximal C-components of the graph
of Figure 1(b) are $\\{X,X_{d},Y,Y_{x},Y_{d}\\}$, $\\{D,D_{x}\\}$,
$\\{Z,Z_{x},Z_{d}\\}$, and $\\{W,W_{x},W_{d}\\}$.

We recall the ID* and IDC* algorithms of Shpitser and Pearl (2007) which are
depicted in Figures 3 and 4 for identifying counterfactual queries and
conditional counterfactual queries, respectively. Both algorithms are sound
and complete (Shpitser and Pearl, 2008, Theorems 26 and 31), meaning that when
they succeed in identifying the query, the expression returned is equal to
$P(\gamma)$ or $P(\gamma|\delta)$, respectively, and when they fail, the query
is not identifiable. We aim to characterize the operation of these algorithms
on an intuitive level and provide line-by-line examples of their operation via
examples.

function ID*($G$, $\gamma$)
---
INPUT: $G$ a causal diagram, $\gamma$ a conjunction of counterfactual events
OUTPUT: an expression for $P(\gamma)$ in terms of $P_{*}$, or FAIL

  1. 1.

if $\gamma=\emptyset$, return 1

  2. 2.

if $(\exists x_{x^{\prime}..}\in\gamma)$, return 0

  3. 3.

if $(\exists x_{x..}\in\gamma)$, return ID*$(G,\gamma\setminus\\{x_{x..}\\})$

  4. 4.

$(G^{\prime},\gamma^{\prime})$ = make-cg$(G,\gamma)$

  5. 5.

if $\gamma^{\prime}$ = INCONSISTENT, return 0

  6. 6.

if $C(G^{\prime})=\\{\mathbf{S}^{1},\ldots,\mathbf{S}^{k}\\}$, return
$\sum_{V(G^{\prime})\setminus\gamma^{\prime}}\prod_{i=1}^{k}$ID*$(G,\mathbf{s}^{i}_{v(G^{\prime})\setminus\mathbf{s}^{i}})$

  7. 7.

if $C(G^{\prime})=\\{\mathbf{S}\\}$, then

    1. 8.

if $(\exists\mathbf{x},\mathbf{x}^{\prime})$ s.t.
$\mathbf{x}\neq\mathbf{x}^{\prime},\mathbf{x}\in\mathrm{sub}(\mathbf{S}),\mathbf{x}^{\prime}\in\mathrm{ev}(\mathbf{S})$,
throw FAIL

    2. 9.

else, let $\mathbf{x}=\cup\,\mathrm{sub}(\mathbf{S})$, return
$P_{\mathbf{x}}(\mathrm{var}(\mathbf{S}))$.

Figure 3: Counterfactual identification algorithm ID* by Shpitser and Pearl
(2007).

We begin by describing the ID* algorithm. On line 1, we check for an empty
conjunction, which by convention has probability 1. Line 2 checks for direct
inconsistencies meaning counterfactual variables that are intervened on but
simultaneously observed to have a different value than the intervention. Such
counterfactuals violate the Axiom of Effectiveness (Pearl, 2009), and if
found, we return probability 0. Line 3 removes tautological counterfactuals
from the conjunction meaning counterfactuals where the variable was observed
to have the value it was forced to take via intervention. Line 4 calls the
make-cg algorithm to construct the counterfactual graph $G^{\prime}$ and the
corresponding conjunction $\gamma^{\prime}$ where some counterfactual
variables may have been relabeled due to equivalence between counterfactual
variables. We leave the details of the make-cg algorithm and the related core
results to Appendix A. In summary, the output $G^{\prime}$ of make-cg is a
refined version of the parallel worlds graph of $G$ and $\gamma$, where each
counterfactual variable is unique. Similarly, if some variables in $\gamma$
were found to be equivalent, then those variables are replaced in
$\gamma^{\prime}$ by their new representatives in $G^{\prime}$. If as a result
of this operation the refined conjunction $\gamma^{\prime}$ is now
inconsistent, we again return probability 0. The next two lines take advantage
of the C-component factorization of the counterfactual graph $G^{\prime}$,
analogously to the ID algorithm. If there is more than one maximal C-component
of $G^{\prime}$, then we proceed to line 6 where the original query is
decomposed into a set of subproblems, each of which we again call ID* for.
Note that the sets $\mathbf{S}^{i}$ are sets of counterfactual variables, but
we may interpret them as counterfactual conjunctions in the subsequent
recursive calls. Similarly, we may interpret $\gamma^{\prime}$ as a set of
counterfactual variables when carrying out the outermost summation over the
possible values of the counterfactual variables in
$V(G^{\prime})\setminus\gamma^{\prime}$. In cases where a set $\mathbf{S}^{i}$
contains counterfactual variables, the intervention
$\textrm{do}(v(G^{\prime})\setminus\mathbf{s}^{i})$ should be understood as
merging of the subscripts, e.g., if $\mathbf{S}^{i}=\\{Y_{x}\\}$ and
$V(G^{\prime})\setminus\mathbf{S}^{i}=\\{Z\\}$, and $Y_{x}$ has the value $y$
in $\gamma^{\prime}$, then
$\mathbf{s}^{i}_{v(G^{\prime})\setminus\mathbf{s}^{i}}=y_{x,z}$.

If there is only one C-component, we enter line 7 that serves as the base
case. There are now only two options. If there is an inconsistent value
assignment on line 8 such that at least one of the values is in the subscript,
then the query is not identifiable, and we fail. If there is no such conflict,
we can take the union of all the subscripts in $\gamma^{\prime}$ and return
their effect on the variables in $\gamma^{\prime}$ on line 9.

function IDC*($G$, $\gamma$, $\delta$)
---
INPUT: $G$ a causal diagram, $\gamma,\delta$ conjunctions of counterfactual
events
OUTPUT: an expression for $P(\gamma|\delta)$ in terms of $P_{*}$, or FAIL, or
UNDEFINED

  1. 1.

if ID*$(G,\delta)=0$, return UNDEFINED

  2. 2.

$(G^{\prime},\gamma^{\prime}\wedge\delta^{\prime})$ = make-
cg$(G,\gamma\wedge\delta)$

  3. 3.

if $\gamma^{\prime}\wedge\delta^{\prime}$ = INCONSISTENT, return $0$

  4. 4.

if $(\exists y_{\mathbf{x}}\in\delta^{\prime})$ s.t.
$(Y_{\mathbf{x}}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\gamma^{\prime})G^{\prime}_{\underline{y_{\mathbf{x}}}}$,
return
IDC*$(G,\gamma^{\prime}_{y_{\mathbf{x}}},\delta^{\prime}\setminus\\{y_{\mathbf{x}}\\})$

  5. 5.

else, let $P^{\prime}$ = ID*$(G,\gamma^{\prime}\wedge\delta^{\prime})$, return
$P^{\prime}/P^{\prime}(\delta)$

Figure 4: Conditional counterfactual identification algorithm IDC* by Shpitser
and Pearl (2007).

In contrast, the IDC* algorithm is simpler, as it leverages the ID* algorithm.
The consistency of the conditioning conjunction $\delta$ is first confirmed on
line 1, and if $\delta$ is found to be inconsistent, then the conditional
probability $P(\gamma|\delta)$ is undefined, and we return. Line 2 applies the
make-cg algorithm to the joint conjunction $\gamma\wedge\delta$ to construct
the corresponding counterfactual graph $G^{\prime}$ and the restructured
version of the conjunction, $\gamma^{\prime}\wedge\delta^{\prime}$. If
$\gamma^{\prime}\wedge\delta^{\prime}$ was found to be inconsistent, we return
probability 0 on line 3. Line 4 takes advantage of conditional independence
relations implied by the counterfactual graph $G^{\prime}$ and the second rule
of do-calculus to add variables as interventions to $\gamma^{\prime}$ by
removing them from $\delta^{\prime}$. If the necessary d-separation holds, we
initiate a recursive call to IDC* again. Finally on line 5, if no more
variables can be removed from $\delta^{\prime}$, we simply apply the ID*
algorithm to the joint conjunction $\gamma^{\prime}\wedge\delta^{\prime}$ and
obtain the identifying functional as a standard conditional probability from
the distribution returned by ID*.

###  3.1 Examples on the identifiability algorithms

We recall the counterfactual conjunction $\gamma=y_{x}\wedge x^{\prime}\wedge
z_{d}\wedge d$ from Section 2 and describe how the ID* algorithm operates when
applied to $P(\gamma)$ in the graph of Figure 1(a), which we will label as $G$
in the context of this example. We start from line 1 and continue to line 2 as
$\gamma$ is not an empty conjunction. On line 2, we note that $\gamma$ does
not contain any inconsistencies, similarly on line 3 we see that $\gamma$ does
not contain any tautological statements. Thus, we reach line 4 and apply the
make-cg algorithm to obtain the counterfactual graph $G^{\prime}$ and the
modified conjunction $\gamma^{\prime}$.

We describe the operation of the make-cg algorithm in this instance. The goal
is to determine which variables in the parallel worlds graph of Figure 1(b)
represent the same variable. We consider all variable pairs in a topological
order of $G$ that originate from the same non-counterfactual variable in $G$.
First, we can conclude that $X$ and $X_{d}$ are the same variable, as they
have the same functional mechanisms and the same parent $U$. By the same
argument, $D$ and $D_{x}$ are the same variable with the common parent
$U_{D}$. The fixed variables $x$ and $d$ cannot be merged with the other
$X$-derived variables and $D$-derived variables, respectively, as their
functional mechanisms are different. Next, we merge $W$ and $W_{d}$ because
their $X$-derived parents ($X$ and $X_{d}$) were found to be the same and they
have the same parent $U_{W}$. However, $W_{X}$ cannot be merged with the other
two $W$-derived variables, because $X$ (and thus $X_{d}$) was observed to
attain the value $x^{\prime}$ in $\gamma$, but $x$ has the value $x$ as fixed
by the intervention. In contrast, we can merge the triplet $Z$, $Z_{x}$ and
$Z_{d}$, because their $D$-derived parents attain the same value, and they
have the same parent $U_{Z}$. The intuition is that because the $U$-variables
are shared between worlds, intervention and observation have the same effect
if the observed values agree with the values fixed by intervention. This is a
consequence of the Axiom of Composition as was considered in the example of
Section 2.1. Finally, we consider the $Y$-derived variables and merge $Y_{x}$
and $Y_{d}$ because their $Z$-derived parents are the same, their $W$-derived
parents are the same, and they have the same parent $U$. The variable $Y_{x}$
cannot be merged with the other two, because its $W$-derived parent $W_{x}$
was not the same variable as $W$ and $W_{d}$.

Consequently, we must choose a name for each merged variable. This choice is
arbitrary and plays no role in the correctness of the algorithm; the
difference is purely notational. In this example, we pick the original name
with the fewest subscripts to represent the merged variable, i.e., $X$
represents the merged pair $X,X_{d}$, $Z$ represents the merged triplet
$Z,Z_{x},Z_{d}$, $W$ represents the merged pair $W,W_{d}$ and finally $Y$
represents the merged pair $Y,Y_{d}$. Note that because the $Z$-derived
variables were all merged but $d$ was not merged with $D$ and $D_{x}$, we
essentially have two $D$-derived parents for the merged $Z$. In such
scenarios, we simply omit the fixed version of the parent variable from the
graph, because this scenario may only arise if the parent variables were found
to have the same value, thus their role in the functional mechanisms of their
children is identical. Lastly, we may restrict our attention to those
counterfactual variables that are ancestral to the query $\gamma$ in this
merged graph, which are $x,W_{x},Y_{x},Z,D,X$ and $U$

Thus, we obtain the counterfactual graph $G^{\prime}$ for $\gamma$ depicted in
Figure 5 using once again the convention that unobserved variables with only
one child are not drawn. As a result of the variable merges, we also update
our original conjunction $\gamma$ with references to the merged variables to
obtain $\gamma^{\prime}=y_{x}\wedge x^{\prime}\wedge z\wedge d$. The new
conjunction $\gamma^{\prime}$ is not inconsistent on line 5, and thus we
continue.

0$\scriptstyle X$0$\scriptstyle D$0$\scriptstyle Z$0$x$0$\scriptstyle
W_{x}$0$\scriptstyle Y_{x}$ Figure 5: Counterfactual graph $G^{\prime}$ for
$y_{x}\wedge x^{\prime}\wedge z_{d}\wedge d$ of the graph of Figure 1(a).

On line 6 we first determine the maximal C-components of the counterfactual
graph $G^{\prime}$ which are $\\{X,Y_{x}\\}$, $\\{Z\\}$, $\\{W_{x}\\}$ and
$\\{D\\}$. By the C-component factorization we have that

$P(y_{x}\wedge x^{\prime}\wedge z\wedge d)=\sum_{w}P(y_{x,z,w,d}\wedge
x^{\prime}_{z,w,d})P(z_{y,x,w,d})P(w_{x,y,z,d})P(d_{y,x,z,w}),$ (1)

which means that we launch four recursive calls to ID* to identify each of the
terms in the right-hand side expression. We will consider the last three terms
first as they result in a similar simple path through the algorithm. For each
of these terms, the counterfactual graph will contain a single non-fixed
vertex ($Z_{y,x,w,d}$, $W_{x,y,z,d}$ and $D_{y,x,z,w}$, respectively). Because
the conjunctions are not empty, there are no inconsistencies or tautologies,
and only a single C-component, we end up on line 7 in each case. None of the
terms contain value assignments that would conflict with the subscript and
thus each term is identified as an interventional distribution on line 9. Note
that when line 7 is reached, redundant subscripts should be removed, i.e.,
those subscript variables that are not ancestors of the counterfactual
variables in $\gamma^{\prime}$ in the counterfactual graph $G^{\prime}$.
Otherwise, a conflict may be found erroneously on line 8. This operation was
not formally included in the algorithm by Shpitser and Pearl (2007), but
nonetheless carried out in a running example by Shpitser and Pearl (2008).
Thus $P(z_{y,x,w,d})=P_{d}(z)$, $P(w_{x,y,z,d})=P_{x}(w)$ and
$P(d_{y,x,z,w})=P(d)$. For the first term $P(y_{x,z,w,d}\wedge
x^{\prime}_{z,w,d})$, the only difference is that the counterfactual graph has
two non-fixed vertices, but the outcome is the same and we end up on line 7
due to the single C-component containing $Y_{x,z,w,d}$ and $X_{z,w,d}$. There
are no conflicts this time either, and we obtain $P(y_{x,z,w,d}\wedge
x^{\prime}_{z,w,d})=P_{w,z}(y,x^{\prime})$. Thus, we obtain the identifying
functional of the counterfactual query:

$P(y_{x}\wedge x^{\prime}\wedge z_{d}\wedge
d)=\sum_{w}P_{w,z}(y,x^{\prime})P_{d}(z)P_{x}(w)P(d).$

Next, we will consider an example that causes a conflict at line 7 resulting
in a non-identifiable counterfactual query. Suppose that we also have an edge
from $X$ to $Y$ in the graph of Figure 1(a) and we wish to identify the same
counterfactual query $P(y_{x}\wedge x^{\prime}\wedge z\wedge d)$ as in the
previous example in this modified graph. The ID* algorithm proceeds similarly
as in the previous example up to line 4 where we obtain a slightly different
counterfactual graph, which is the graph of Figure 5, but with the
corresponding extra edge from $X$ to $Y_{x}$. Thus, the algorithm proceeds
similarly to line 6, where the C-component factorization is the same as (1).
The last three terms are still identifiable, but this time the first term
$P(y_{x,z,w,d}\wedge x^{\prime}_{z,w,d})$ is problematic. On line 7 after
removing redundant interventions, the term takes the form $P(y_{x,z,w}\wedge
x^{\prime}_{z,w})$ which now contains a conflict, because $x$ appears in the
subscript but $x^{\prime}$ is observed at the same time, resulting in non-
identification on line 8.

We return to the example presented in Section 2.1 and apply the IDC* algorithm
to identify the counterfactual query $P(y_{x}|z_{x}\wedge x^{\prime})$ in the
graph of Figure 2, which we will again refer to as $G$ in the context of this
example. We trace the application of IDC*$(G,y_{x},z_{x}\wedge x^{\prime})$.
On line 1, the ID* algorithm is applied to $z_{x}\wedge x^{\prime}$, which is
not identifiable, but also not inconsistent. Continuing to line 2, we apply
the make-cg algorithm to construct the counterfactual graph $G^{\prime}$,
which is shown in Figure 6(a).

0$\scriptstyle X$0$x$0$\scriptstyle Z$0$\scriptstyle Z_{x}$0$\scriptstyle
Y$0$\scriptstyle Y_{x}$0$\scriptstyle U$0$\scriptstyle U_{Y}$ (a) Parallel
worlds graph for $y_{x}\wedge z_{x}\wedge x^{\prime}$ (the counterfactual
graph).

0$\scriptstyle X$0$x$0$\scriptstyle Z$0$z$0$\scriptstyle Y$0$\scriptstyle
Y_{x,z}$0$\scriptstyle U$0$\scriptstyle U_{Y}$ (b) Parallel worlds graph for
$y_{x,z}\wedge x^{\prime}$ (the counterfactual graph).

Figure 6: Counterfactual graphs used during the derivation of
$P(y_{x}|z_{x}\wedge x^{\prime})$.

Because $X$ was observed to have the value $x^{\prime}$, but the intervention
for $Z$ and $Y$ has the value $x$, we cannot merge $X$ and $x$. Similarly, the
$X$-parent of $Z$ in both worlds has a different value, meaning that $Z$ and
$Z_{x}$ cannot be merged either. Finally, through the same reasoning, $Y$ and
$Y_{x}$ will remain unmerged due to the difference in the $Z$-parent. Thus,
the parallel worlds graph is the counterfactual graph $G^{\prime}$ in this
instance. This also means that $\gamma^{\prime}=\gamma$ and
$\delta^{\prime}=\delta$ in the output of make-cg.

On line 3, we check for inconsistencies in $y_{x}\wedge z_{x}\wedge
x^{\prime}$, but there are none. Next on line 4, we check whether either of
the two variables in $\delta^{\prime}$ are d-separated from $\gamma^{\prime}$
when outgoing edges of that variable have been removed. We can see that $X$ is
not d-separated from $Y_{x}$, because the path $X\leftarrow U_{1}\rightarrow
Z_{x}\rightarrow Y_{x}$ is open in $G^{\prime}_{\underline{X}}$. However,
$Z_{x}$ is d-separated from $Y_{x}$ in $G^{\prime}_{\underline{Z_{x}}}$ (note
that $x$ is fixed by intervention, and thus the path $Z_{x}\leftarrow
x\rightarrow Y_{x}$ is not an open backdoor path). Thus, line 4 adds an
intervention on $Z$ to $Y_{x}$ because $Y_{x}$ is a descendant of $Z_{x}$ in
$G^{\prime}$, and removes $Z_{x}$ from $\delta^{\prime}$, and we call
IDC*$(G^{\prime},y_{x,z},x^{\prime})$.

We now trace this new recursive call. Once again on line 1, ID* is not able to
identify the effect, but is also not inconsistent. Next, we construct a new
counterfactual graph $G^{\prime\prime}$ for $y_{x,z}\wedge x^{\prime}$ as
depicted in Figure 6(b). Using similar reasoning as before, the make-cg
algorithm is not able to merge any nodes this time either and thus the
parallel worlds graph is the counterfactual graph. Again, this means that
$\gamma^{\prime\prime}=\gamma^{\prime}$ and
$\delta^{\prime\prime}=\delta^{\prime}$ in the output of make-cg. Line 3
checks again for inconsistencies in $y_{x,z}\wedge x^{\prime}$, but there are
none. Thus we arrive again on line 4, but this time $X$ is d-separated from
$Y_{x,z}$ in $G^{\prime\prime}_{\underline{X}}$. Now, $Y_{x,z}$ is not a
descendant of $X$ in $G^{\prime\prime}$ so no new intervention is added to
$Y_{x,z}$, and $x^{\prime}$ is removed from $\delta^{\prime\prime}$. Because
the conditioning $\delta$-argument of the next IDC* call is now empty, we can
call ID* directly as ID*$(G,y_{x,z})$, but $P(y_{x,z})$ is no longer a
counterfactual quantity, but an interventional distribution and thus directly
identifiable from $P_{*}$ as $P_{x,z}(y)$.

We note the difference compared to the manual identification strategy we used
in Section 2.1 to obtain identifiability. Instead of using axioms of
counterfactuals or independence restrictions explicitly, the ID* and IDC*
algorithms take full advantage of the counterfactual graph and the conditional
independence relations between the counterfactual variables implied by it.

##  4 The cfid package

The cfid package is available from CRAN at
https://cran.r-project.org/package=cfid and can be obtained in R using the
following commands: R> install.packages("cfid") R> library("cfid") Development
of cfid takes place on GitHub https://github.com/santikka/cfid.

The main contributions of the cfid package are the implementations of the ID*
and IDC* algorithms. The package also provides reimplementations of the ID and
IDC algorithms for interventional distributions from the causaleffect package,
but without relying on the igraph (Csardi and Nepusz, 2006) package. In fact,
cfid has no mandatory package dependencies or installation requirements. The
cfid package provides its own text-based interface for defining graphs, which
closely follows the syntax of the dagitty package, and also supports other
external graph formats directly. Installation of the igraph and dagitty
packages is optional and required only if the user wishes to import or export
graphs using the aforementioned packages.

The inclusion of the identifiability algorithms for interventional
distributions enables a full identification pipeline. First, we determine the
identifiability of a counterfactual query from the set of all interventional
distributions, and then proceed to identify each interventional distribution
that appears in the identifying functional of the counterfactual from the
joint observed probability distribution of the causal model. The level of
attempted identification can be specified by the user.

###  4.1 Defining causal diagrams

Causal diagrams (i.e., DAGs) in cfid are constructed via the function dag
dag(x, u = character(0L)) where x is a single character string in a syntax
analogous to the DOT language for GraphViz (and the dagitty package), and u is
an optional character vector of variable names that should be considered
unobserved in the graph. Internally, a semi-Markovian representation is always
used for DAGs where each latent variable has at most two children, which is
obtained from the input via the latent projection (Verma and Pearl, 1990).

As an example, the graph of Figure 2 can be constructed as follows: R> g <\-
dag("X -> Z -> Y; X -> Y; X <-> Z") Above, individual statements are separated
by a semicolon for additional clarity, but this is optional, and a space would
suffice. More generally, the input of dag consists of statements of the form
$n_{1}e_{1}n_{2}e_{2}\cdots e_{k}n_{k}$ where each $e_{i}$ symbol must be a
supported edge type, i.e., ->, <\- or <->, and each $n_{i}$ symbol must
correspond to single node such as X or a subgraph such as {X, Y, Z} or {X ->
Y}. Subgraphs are enclosed within curly braces, and they follow the same
syntax as x. Subgraphs can also be nested arbitrarily. An edge of the form X
-> {…} means that there is an edge from X to all vertices in the subgraph, and
the interpretation for <\- and <-> is analogous. Individual statements in the
graph definition can be separated by a semicolon, a space, or a new line.
Commas can be used within subgraphs to distinguish vertices, but a space is
sufficient.

The same DAG can often be defined in many ways. For example, we could also
define the graph of Figure 2 using a subgraph construct as follows: R> g <\-
dag("X -> Z, Y; Z -> Y; X <-> Z") We could also combine the outgoing edge of Z
and the bidirected edge into a single statement: R> g <\- dag("X -> Z, Y; X
<-> Z -> Y;") The edge from Z to Y could be defined in the subgraph as well:
R> g <\- dag("Z <-> X -> Z -> Y") The output of dag is an object of class
"dag" which is a square adjacency matrix of the graph, with additional
attributes for the vertex labels and latent variables and a print method.
Graph definitions that imply cycles or self-loops will raise an error.
Examples of more complicated graph constructs can be found from the cfid
package documentation for the dag function. Graphs using supported external
formats can be converted to dag objects via the function import_graph.
Conversely, dag objects can be exported in supported external formats using
the function export_graph.

###  4.2 Defining counterfactual variables and conjunctions

Counterfactual variables are defined via the function counterfactual_variable
or its shorthand alias cf counterfactual_variable(var, obs = integer(0L), sub
= integer(0L)) cf(var, obs = integer(0L), sub = integer(0L)) The first
argument var is a single character string naming the variable, e.g., "Y". The
second argument obs describes the value assignment as a single integer. The
value of this argument does not describe the actual value taken by the
variable, but simply the assignment level, meaning that obs = 1 is a different
value assignment than obs = 0, but the actual values that the counterfactual
variable takes need not necessarily be 1 and 0. The idea is similar to the
internal type of factors in R. Finally, sub defines the set of interventions
as a named integer vector, where the actual values correspond to the
intervention levels, and not actual values, analogous to obs. The output of cf
is an object of class "counterfactual_variable".

As an example, the counterfactual variables in $\gamma=y_{x}\wedge
x^{\prime}\wedge z_{d}\wedge d$ can be defined as follows: R> v1 <\- cf(var =
"Y", obs = 0L, sub = c(X = 0L)) R> v2 <\- cf(var = "X", obs = 1L) R> v3 <\-
cf(var = "Z", obs = 0L, sub = c(D = 0L)) R> v4 <\- cf(var = "D", obs = 0L) R>
list(v1, v2, v3, v4) [[1]] y_x

[[2]] x’

[[3]] z_d

[[4]] d The print method for counterfactual_variable objects mimics the
notation used in this paper in LaTeX syntax.

Individual counterfactual_variable objects can be combined into a
counterfactual conjunction via the function counterfactual_conjunction or its
shorthand alias conj. This function takes arbitrarily many
counterfactual_variable objects as input. The output of conj is an object of
class "counterfactual_conjunction". R> c1 <\- conj(v1, v2, v3, v4) R> c1 y_x /
x’ / z_d / d Alternatively, the ‘+‘ operator can be used to build conjunctions
from counterfactual variables or conjunctions. R> c2 <\- v1 + v2 R> c3 <\- v3
+ v4 R> c2 R> c3 R> c2 + c3 y_x / x’ z_d / d y_x / x’ / z_d / d The subset
operator ‘[‘ is supported for counterfactual conjunctions R> c1[c(1, 3)] y_x /
z_d Just as the cf function, the print method for counterfactual_conjunction
objects mimics the formal notation of using the $\wedge$ symbol to separate
individual statements, but this symbol can also be changed by the user.

###  4.3 Identifying counterfactual queries

Identification of counterfactual queries is carried out by the function
identifiable identifiable(g, gamma, delta = NULL, data = "interventions")
where g is a causal diagram defined by the function dag, gamma is the
conjunction $\gamma$ as a counterfactual_conjunction object describing the
counterfactual query $P(\gamma)$ to be identified, delta is an optional
argument also of class counterfactual_conjunction that should be provided if
identification of a conditional counterfactual $P(\gamma|\delta)$ is desired
instead. Finally, data defines the available probability distributions for
identification. The default value "interventions" means that identification is
carried out to the intervention level, i.e., by using only the set of all
interventional distributions $P_{*}$. The alternatives are "observations",
where only the joint observed probability distribution $P(\mathbf{v})$ is
available, and "both" where both $P_{*}$ and $P(\mathbf{v})$ are available,
and identification in terms of $P(\mathbf{v})$ is prioritized.

We reassess the identifiability examples of Section 3.1 using the cfid
package. The conjunction of the query $\gamma$ for the first two examples has
already been defined as c1 in the previous section. We define the graphs for
the identifiable case in Figure 1(a) and the non-identifiable case with the
additional edge from $X$ to $Y$: R> g1 <\- dag("Y <-> X -> W -> Y <\- Z <\-
D") R> g2 <\- dag("Y <-> X -> W -> Y <\- Z <\- D; X -> Y") R> out1 <\-
identifiable(g1, c1) R> out2 <\- identifiable(g2, c1) R> out1 R> out2 The
query P(y_x / x’/ z_d / d) is identifiable from P_*. Formula: ∑_w
P_w,z(y,x’)P_x(w)P_d(z)P(d)

The query P(y_x / x’/ z_d / d) is not identifiable from P_*. The identifiable
function returns an object of class "query", whose print method provides a
summary of the identification result. Objects of this class are lists with the
following elements:

id

A logical value that is TRUE is the counterfactual query is identifiable and
FALSE otherwise

formula

An object of class "functional" representing the identifying functional. The
format method for functional objects provides the formula of the
counterfactual query in LaTeX syntax when the query is identifiable.
Otherwise, formula is NULL.

undefined

A logical value that is TRUE is a conditional counterfactual query is found to
be undefined.

query

The original query as a counterfactual_conjunction object.

data The data argument passed to identifiable.

By default, the notation of Shpitser and Pearl (2007) is used for
interventional distributions, where interventions are denoted using the
subscript, e.g. $P_{x}(y)$. If desired, the notation can be swapped to Pearl’s
notation with the explicit do-operator denoting interventions, e.g.,
$P(y|\textrm{do}(x))$. This can be accomplished via the use_do argument of the
format method for functional objects (passed here via the print method): R>
print(out1[["formula"]], use_do = TRUE) ∑_w
P(y,x’|do(w,z))P(w|do(x))P(z|do(d))P(d) For the third example of Section 3.1,
we have already defined the counterfactual variable $Y_{x}$ of the query as v1
and the observation $X=x^{\prime}$ in the condition as v2. We still need to
define the graph of Figure 2 and the other conditioning variable $Z_{x}$: R>
g3 <\- dag("Z <-> X -> Z -> Y") R> v5 <\- cf("Z", 0, c(X = 0)) R>
identifiable(g3, v1, v5 + v2) The query P(y_x|z_x / x’) is identifiable from
P_*. Formula: P_x,z(y) Recall from Section 2.1, that this interventional
distribution can be further identified, which can be accomplished by setting
the data argument to "observations" in identifiable (or to "both" in this
case): R> identifiable(g3, v1, v5 + v2, data = "observations") The query
P(y_x|z_x / x’) is identifiable from P(v). Formula: P(y|x,z)

##  5 Summary

The cfid package provides an easy-to-use interface to identifiability analysis
of counterfactual queries. The causal diagram of the causal model can be
specified by the user via an intuitive interface, and a variety commonly used
external graph formats are supported. The results from the identifiability
algorithms are wrapped neatly in a LaTeX syntax to be readily used in
publications or reports. This tutorial demonstrates the features of the
package and provides insight into the core algorithms it implements.

## Acknowledgments

This work was supported by Academy of Finland grant number 331817.

## References

  * Avin _et al._ (2005) Avin C, Shpitser I, Pearl J (2005).  “Identifiability of Path-Specific Effects.”  In _Proceedings of International Joint Conference on Artificial Intelligence_ , volume 19, pp. 357–363.
  * Balke and Pearl (1994a) Balke A, Pearl J (1994a).  “Counterfactual Probabilities: Computational Methods, Bounds and Applications.”  In _Proceedings of the 10th Conference on Uncertainty in Artificial Intelligence_ , pp. 46–54.
  * Balke and Pearl (1994b) Balke A, Pearl J (1994b).  “Probabilistic Evaluation of Counterfactual Queries.”  In _Proceedings of the 12th AAAI National Conference on Artificial Intelligence_ , pp. 230–237.
  * Bareinboim and Pearl (2012) Bareinboim E, Pearl J (2012).  “Causal Inference by Surrogate Experiments: $z$-Identifiability.”  In _Proceedings of the 28th Conference on Uncertainty in Artificial Intelligence_ , pp. 113–120.
  * Bhattacharya _et al._ (2020) Bhattacharya R, Nabi R, Shpitser I (2020).  “Semiparametric Inference For Causal Effects In Graphical Models With Hidden Variables.”  10.48550/ARXIV.2003.12659.  URL https://arxiv.org/abs/2003.12659.
  * Chen _et al._ (2020) Chen M, Chernozhukov V, Fernandez-Val I, Melly B (2020).  _Counterfactual: Estimation and Inference Methods for Counterfactual Analysis_.  R package version 1.2, URL https://CRAN.R-project.org/package=Counterfactual.
  * Csardi and Nepusz (2006) Csardi G, Nepusz T (2006).  “The igraph Software Package for Complex Network Research.”  _InterJournal_ , Complex Systems, 1695.  URL https://igraph.org.
  * Halpern (1998) Halpern JY (1998).  “Axiomatizing Causal Reasoning.”  In _Proceedings of the 14th Conference on Uncertainty in Artificial Intelligence_ , pp. 202–210.
  * Holland (1986) Holland PW (1986).  “Statistics and Causal Inference.”  _Journal of the American Statistical Association_ , 81(396), 945–960.  10.1080/01621459.1986.10478354.
  * Huang and Valtorta (2006) Huang Y, Valtorta M (2006).  “Pearl’s Calculus of Intervention is Complete.”  In _Proceedings of the 22nd Conference on Uncertainty in Artificial Intelligence_ , pp. 217–224. AUAI Press.
  * Karvanen (2022) Karvanen J (2022).  _R6causal: R6 Class for Structural Causal Models_.  R package version 0.6.1.
  * Kivva _et al._ (2022) Kivva Y, Mokhtarian E, Etesami J, Kiyavash N (2022).  “Revisiting the General Identifiability Problem.”  In _Proceedings of the 38th Conference on Uncertainty in Artificial Intelligence_ , volume 180, pp. 1022–1030. PMLR.
  * Kusner _et al._ (2017) Kusner MJ, Loftus J, Russell C, Silva R (2017).  “Counterfactual Fairness.”  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , pp. 4069–4079.
  * Lee and Shpitser (2020) Lee JJR, Shpitser I (2020).  “Identification Methods With Arbitrary Interventional Distributions as Inputs.”  10.48550/ARXIV.2004.01157.  URL https://arxiv.org/abs/2004.01157.
  * Lee _et al._ (2019) Lee S, Correa JD, Bareinboim E (2019).  “General Identifiability with Arbitrary Surrogate Experiments.”  In _Proceedings of the 35th Conference on Uncertainty in Artificial Intelligence_ , volume 115, pp. 389–398. PMLR.
  * Nabi _et al._ (2020) Nabi R, Bhattacharya R, Shpitser I (2020).  “Full Law Identification in Graphical Models of Missing Data: Completeness Results.”  In _Proceedings of the 37th International Conference on Machine Learning_ , volume 119, pp. 7153–7163.  10.5555/3524938.3525601.
  * Pearl (1995) Pearl J (1995).  “Causal Diagrams for Empirical Research.”  _Biometrika_ , pp. 669–710.  10.1093/biomet/82.4.669.
  * Pearl (2009) Pearl J (2009).  _Causality: Models, Reasoning and Inference_.  2nd edition. Cambridge University Press.
  * R Core Team (2022) R Core Team (2022).  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing.  URL https://www.R-project.org/.
  * Sharma _et al._ (2019) Sharma A, Kiciman E, _et al._ (2019).  “DoWhy: A Python Package for Causal Inference.”  URL https://github.com/microsoft/dowhy.
  * Shpitser and Pearl (2006a) Shpitser I, Pearl J (2006a).  “Identification of Conditional Interventional Distributions.”  In _Proceedings of the 22nd Conference on Uncertainty in Artificial Intelligence_ , pp. 437–444. AUAI Press.
  * Shpitser and Pearl (2006b) Shpitser I, Pearl J (2006b).  “Identification of Joint Interventional Distributions in Recursive Semi-Markovian Causal Models.”  In _Proceedings of the 21st National Conference on Artificial Intelligence - Volume 2_ , pp. 1219–1226. AAAI Press.
  * Shpitser and Pearl (2007) Shpitser I, Pearl J (2007).  “What Counterfactuals Can Be Tested.”  In _Proceedings of the 23rd Conference on Uncertainty in Artificial Intelligence_ , pp. 352–359. AUAI Press.
  * Shpitser and Pearl (2008) Shpitser I, Pearl J (2008).  “Complete Identification Methods for the Causal Hierarchy.”  _Journal of Machine Learning Research_ , 9(64), 1941–1979.
  * Stoll _et al._ (2020) Stoll H, King G, Zeng L, Gandrud C, Sabath B (2020).  _WhatIf: Software for Evaluating Counterfactuals_.  R package version 1.5-10, URL https://CRAN.R-project.org/package=WhatIf.
  * Textor _et al._ (2017) Textor J, van der Zander B, Gilthorpe MS, Liśkiewicz M, Ellison GT (2017).  “Robust Causal Inference Using Directed Acyclic Graphs: The R package dagitty.”  _International Journal of Epidemiology_ , 45(6), 1887–1894.  10.1093/ije/dyw341.
  * Tian and Pearl (2002) Tian J, Pearl J (2002).  “A General Identification Condition for Causal Effects.”  In _Proceedings of the 19th AAAI National Conference on Artificial Intelligence_ , pp. 567–573.
  * Tikka _et al._ (2021) Tikka S, Hyttinen A, Karvanen J (2021).  “Causal Effect Identification from Multiple Incomplete Data Sources: A General Search-Based Approach.”  _Journal of Statistical Software_ , 99(5), 1–40.  10.18637/jss.v099.i05.
  * Tikka and Karvanen (2017) Tikka S, Karvanen J (2017).  “Identifying Causal Effects with the R Package causaleffect.”  _Journal of Statistical Software_ , 76(12), 1–30.  10.18637/jss.v076.i12.  URL https://www.jstatsoft.org/index.php/jss/article/view/v076i12.
  * Tikka and Karvanen (2019) Tikka S, Karvanen J (2019).  “Surrogate Outcomes and Transportability.”  _International Journal of Approximate Reasoning_ , 108, 21–37.
  * Van Rossum and Drake (2009) Van Rossum G, Drake FL (2009).  _Python 3 Reference Manual_.  CreateSpace.
  * Verma and Pearl (1990) Verma TS, Pearl J (1990).  “Equivalence and Synthesis of Causal Models.”  In _Proceedings of the 6th Conference on Uncertainty in Artificial Intelligence_ , pp. 255–270.
  * Zhang and Bareinboim (2018) Zhang J, Bareinboim E (2018).  “Fairness in Decision-Making — The Causal Explanation Formula.”  In _Proceedings of the 32nd AAAI Conference on Artificial Intelligence_ , pp. 2037–2045.

##  Appendix A Counterfactual graphs

We restate here the make-cg algorithm and the associated Lemmas that are used
to construct counterfactual graphs from parallel worlds graphs. Lemma 1 is
used to characterize conditions where two counterfactual variables in fact
represent the same random variable. Lemma 2 shows that under the conditions of
Lemma 1 such variables can be merged into a single random variable. For more
details, see (Shpitser and Pearl, 2008).

######  Lemma 1 (Lemma 24 of Shpitser and Pearl (2008)).

Let $M$ be a model inducing $G$ containing variables $\alpha,\beta$ with the
following properties:

  * •

$\alpha$ and $\beta$ have the same domain of values.

  * •

There is a bijection $f$ from $Pa(\alpha)$ to $Pa(\beta)$ such that a parent
$\gamma$ and $f(\gamma)$ have the same domain of values.

  * •

The functional mechanisms of $\alpha$ and $\beta$ are the same (except
whenever the function for $\alpha$ uses the parent $\gamma$, the corresponding
function for $\beta$ uses $f(\gamma)$.

Assume as observable variable set $\mathbf{Z}$ was observed to attain values
$\mathbf{z}$ in $M_{\mathbf{x}}$, the submodel obtained from $M$ by forcing
another observable variable set $\mathbf{X}$ to attain values $\mathbf{x}$.
Assume further that for each $\gamma\in Pa(\alpha)$, either
$f(\gamma)=\gamma$, or $\gamma$ and $f(\gamma)$ attain the same values
(whether by observation or intervention). Then $\alpha$ and $\beta$ are the
same random variable in $M_{\mathbf{x}}$ with observations $\mathbf{z}$.

######  Lemma 2 (Lemma 25 of Shpitser and Pearl (2008)).

Let $M_{\mathbf{x}}$ be a submodel derived from $M$ with set $\mathbf{Z}$
observed to attain values $\mathbf{z}$, such that Lemma 1 holds for
$\alpha,\beta$. Let $M^{\prime}$ be a causal model obtained from $M$ by
merging $\alpha,\beta$ into a new node $\omega$, which inherits all parents
and the functional mechanism of $\alpha$. All children of $\alpha,\beta$ in
$M^{\prime}$ become children of $\omega$. Then $M_{\mathbf{x}}$,
$M^{\prime}_{\mathbf{x}}$ agree on any distribution consistent with
$\mathbf{z}$ being observed.

The previous two Lemmas are leveraged in the make-cg algorithm as shown in
Figure 7. In the implementation of this algorithm, equivalence classes of
worlds for each variable are initialized such that each world is the only
member of its equivalence class at the beginning. During the iteration of step
3, we update these equivalence classes such that when two instances of the
same original variable are found to be the same variable, we combine the
equivalence classes of the worlds that the two instances $\alpha$ and $\beta$
of this variable belong to. When applying Lemma 1, it is not necessary to
check the values of the unobserved parents for equality, because unobserved
parents are shared between worlds, and thus their values will always be equal
when two instances of the same variable are compared.

In the implementation, all graph modifications are carried out at once,
because it is not necessary to modify the graph when determining which
variable pairs are equivalent. In other words, in the implementation of step
3, we only iterate step 3.3 while dynamically updating a list of variable
merges which are all carried out at once before step 4 takes place (i.e., all
modifications specified at steps 3.1 and 3.2). This dynamic approach also
allows us to skip some redundant comparisons. For example, say that variables
$A$ and $B$ have been found to be equivalent, now we only need to compare
either $A$ or $B$ to a third variable $C$.

function make-cg($G$, $\gamma$)
---
INPUT: $G$ a causal diagram, $\gamma$ a conjunction of counterfactual events
OUTPUT: A counterfactual graph $G_{\gamma}$, and either a set of events
$\gamma^{\prime}$ s.t. $P(\gamma^{\prime})=P(\gamma)$ or
INCONSISTENT

  1. 1.

Construct a submodel graph $G_{\mathbf{X}_{i}}$ for each action
$do(\mathbf{x}_{i})$ mentioned in $\gamma$. Construct the parallel worlds
graph $G^{\prime}$ by having all such submodel graphs share their
corresponding $U$ nodes.

  2. 2.

Let $\pi$ be a topological ordering of nodes in $G^{\prime}$, let
$\gamma^{\prime}:=\gamma$.

  3. 3.

Apply Lemmas 1 and 2, in order $\pi$, to each observable node pair
$\alpha,\beta$ derived from the same variable in $G$. For each $\alpha,\beta$
that are the same do:

    1. 3.1

Let $G^{\prime}$ be modified as specified in Lemma 2.

    2. 3.2

Modify $\gamma^{\prime}$ by renaming all occurrences of $\beta$ to $\alpha$.

    3. 3.3

If $\mathrm{val}(\alpha)\neq\mathrm{val}(\beta)$, return ($G^{\prime}$,
INCONSISTENT).

  4. 4.

Return $(G^{\prime}_{An(\gamma^{\prime})},\gamma^{\prime})$, where
$An(\gamma^{\prime})$ is the set of nodes in $G^{\prime}$ ancestral to nodes
corresponding to variables mentioned in $\gamma^{\prime}$.

Figure 7: An algorithm for constructing counterfactual graphs.