Wang et al. RRd stars as robust distance and metallicity indicators

17 2023 10.1017/xxxxx

Proceedings IAU Symposium

# Double-mode RR Lyrae star — robust distance and metallicity indicators

Shu Wang1,3 Xiaodian Chen1,2,3 Jianxing Zhang1,2 & Licai Deng1,2,3 1CAS Key
Laboratory of Optical Astronomy, National Astronomical Observatories, Chinese
Academy of Sciences, Beijing 100101, China
email<EMAIL_ADDRESS>
2School of Astronomy and Space Science, University of the Chinese Academy of
Sciences, Beijing, 100049, China
3Department of Astronomy, China West Normal University, Nanchong, 637009,
China

###### Abstract

RR Lyrae (RR Lyr) stars are a well-known and useful distance indicator for old
stellar populations such as globular clusters and dwarf galaxies. Fundamental-
mode RR Lyr (RRab) stars are commonly used to measure distances, and the
accuracy of the determined distance is strongly constrained by metallicity.
Here, we investigate the metallicity dependence in the period–luminosity (PL)
relation of double-mode RR Lyr (RRd) stars. We find and establish a linear
relation between metallicity and period or period ratio for RRd stars. This
relation can predict the metallicity as accurately as the low-resolution
spectra. Based on this relation, we establish a metallicity-independent PL
relation for RRd stars. Combining the distance of the Large Magellanic Cloud
and Gaia parallaxes, we calibrate the zero point of the derived PL relation to
an error of 0.022 mag. Using RRd stars, we measure the distances of globular
clusters and dwarf galaxies with an accuracy of 2-3% and 1-2%, respectively.
In the future, RRd stars could anchor galaxy distances to an accuracy of 1.0%
and become an independent distance ladder in the Local Group.

###### keywords:

RR Lyrae variable stars, RRd variable stars, distance indicators, metallicity,
dwarf galaxies

††editors: R. de Grijs, P. A. Whitelock & M. Catelan, eds.

##  1 Introduction

RR Lyrae (RR Lyr) stars are low-mass stars that have evolved to the horizontal
branch stage and begin to pulsate as they move into the instability strip. The
pulsation periods of RR Lyr stars are about 0.2 to 1.0 days. According to the
pulsation mode, RR Lyr stars are classified into RRab (fundamental mode), RRc
(first-overtone mode) and RRd (double or multiple mode) stars. In the
Hertzsprung–Russell diagram, RRd stars are located in the intersection region
of the instability strips of RRab and RRc stars.

Most RRd stars are clustered in a sequence on the Petersen diagram (Petersen,
1973) and are called classical RRd stars. There are also anomalous and
peculiar RRd stars. The anomalous or peculiar RRd stars are located above or
below the sequence of classical RRd stars on the Petersen diagram (Nemec &
Moskalik, 2021). The dominant mode of the anomalous or peculiar RRd stars is
the fundamental mode, while the dominant mode of classical RRd stars is the
first-overtone mode. In addition, anomalous or peculiar RRd stars usually have
long-term amplitude modulation. In this paper, the default RRd star is the
classical RRd star.

RR Lyr stars were searched for by a large number of time domain survey,
including the Optical Gravitational Lensing Experiment (OGLE, Soszyński et
al., 2009; Soszynśki et al., 2016, 2019), the Gaia (Clementini et al., 2022),
the Zwicky Transient Facility (ZTF, Chen et al., 2020), the All-Sky Automated
Survey for Supernovae (ASASSN, Jayasinghe et al., 2018), the Asteroid
Terrestrial-impact Last Alert System (ATLAS, Heinze et al., 2018), the Wide-
field Infrared Survey Explorer (WISE, Chen et al., 2018), and the
PanSTARRS1(Sesar et al., 2017).

RR Lyr stars are good distance tracers for the old population. They satisfy
the metallicity–luminosity (LZ) relation in the $V$ band and the
period–metallicity–luminosity (PLZ) relation in infrared bands. Since the two
periods of RRd stars also satisfy the PLZ relation of RRab or RRc,
respectively, RRd stars also have LZ and PLZ relations. Usually, RRd stars
were not adopted to measure distances because the determined periods of RRd
stars are not as accurate as those of RRab stars.

For RRd stars, both theory (Popielski et al., 2000; Marconi et al., 2015) and
observations (Catelan & Smith, 2015; Soszynśki et al., 2019) suggest that
compared to metal-poor RRd stars, metal-rich RRd stars have shorter
fundamental periods and smaller period ratios. There seems to be a good linear
relation between the period ratio and metallicity (Braga et al., 2022). In
this work, we investigate the metallicity dependence in the period–luminosity
(PL) relation of RRd stars and establish their PL relation (Chen et al.,
2023).

##  2 Period–Metallicity Relation

We collect RRd stars mainly based on Gaia and OGLE databases. From Gaia DR3,
we obtain 1021 Galactic RRd stars. From OGLE, we collect 2083 and 674 RRd
stars belonging to the Large Magellanic Cloud (LMC) and the Small Magellanic
Cloud (SMC), respectively. Metallicity information for these stars is taken
from low-resolution spectra of the Sloan Digital Sky Survey (SDSS, Eisenstein
et al., 2011) and the Large Sky Area Multi-Object Fiber Spectroscopic
Telescope (LAMOST, Cui et al., 2012). After cross-matching, we obtain 100 RRd
stars with SDSS or LAMOST metallicities. In addition, we supplement another
203 RRd stars with SDSS or LAMOST metallicity based on the ZTF data. The final
sample contains 207 and 96 RRd stars with metallicity measurements from SDSS
and LAMOST, respectively.

$\displaystyle{\rm[Fe/H]_{SDSS}}=$ $\displaystyle-(173\pm 38)(\log\frac{P_{\rm
1O}}{P_{\rm F}}-\log 0.745)$ (1) $\displaystyle-(6.62\pm 0.88)(\log P_{\rm
1O}-\log 0.37)-(1.76\pm 0.01),$ $\displaystyle{\rm[Fe/H]_{LAMOST}}=$
$\displaystyle-(184\pm 70)(\log\frac{P_{\rm 1O}}{P_{\rm F}}-\log 0.745)$
$\displaystyle-(5.69\pm 1.50)(\log P_{\rm 1O}-\log 0.37)-(1.68\pm 0.03),$

Based on the selected RRd stars, we establish the period–period
ratio–metallicity relation. The results are shown in equation 1. For the SDSS
metallicity, the scatter $\sigma$ is 0.16 dex and the correlation coefficient
$R^{2}$ is 0.74. For the LAMOST metallicity, the scatter $\sigma$ is 0.21 dex
and the correlation coefficient $R^{2}$ is 0.60. The scatter of this relation
decreases when a tighter criterion on metallicity error is adopted. For
example, when we restrict the metallicity error to less than 0.04, with the
rest of 56 RRd, the scatter of this relation will decrease to 0.13 dex. By
comparing the scatter of the period–period ratio–metallicity relation with the
external error of metallicity, we find that the metallicity estimated from the
period and the period ratio of RRd stars can be as accurate as the low-
resolution spectra.

With metallicities from high-resolution spectra, the scatter of the
period–period ratio–metallicity relation could be greatly reduced. For
example, for Cepheids, the scatter of the relation is only 0.033 dex (Kovtyukh
et al., 2016). In the future, with metallicities from high-resolution spectra,
the period–period ratio–metallicity relation of RRd stars could also have a
small scatter. Kovtyukh et al. (2016) found that for double–mode Cepheids,
their metallicity depend on period ratio more than the period. Based on our
sample, we find that RRd metallicity is more dependent on the period than the
period ratio.

Figure 1: The Period–metallicity diagram of RRd stars with metallicities from
SDSS (green dots) and LAMOST (black pentagrams). The green dashed line and the
black solid line are linear fit lines. Red squares are the grid points of the
theoretical model.

Figure 1 shows the distribution of RRd stars in the period–metallicity (PZ)
diagram. The green dots and black pentagrams are RRd stars with metallicities
from SDSS and LAMOST, respectively. Red squares are the grid points of the
theoretical model (Marconi et al., 2015). There is an intuitive linear
relationship between the metallicity and the period. The green dashed line and
the black solid line are linear fit lines. The determined PZ relations for RRd
stars with SDSS metallicities are listed in equation 2.

$\begin{split}{\rm[Fe/H]_{SDSS}}=(-10.33\pm 0.47)\log P_{\rm F}+(-4.87\pm
0.14),\sigma=0.17\,{\rm dex},R^{2}=0.71,\\\ {\rm[Fe/H]_{SDSS}}=(-10.15\pm
0.46)\log P_{\rm 1O}+(-6.11\pm 0.19),\sigma=0.17\,{\rm
dex},R^{2}=0.71\end{split}$ (2)

##  3 Period–Luminosity Relation

Combining the PLZ relation with the period–period ratio–metallicity relation,
the period–period ratio–luminosity relation is derived. We determined the
period–period ratio–luminosity relation based on OGLE LMC RRd stars. We find
that there is an approximate linear relation between the period and the period
ratio, with a correlation coefficient of 0.7 and a dispersion of 0.0003.
Therefore, we try to simplify the period–period ratio–luminosity relation to
the period–luminosity (PL) relation. We find that the absolute magnitude
difference predicted by the PL relation and the period–period ratio–luminosity
relation is small. The scatter in magnitude difference is 0.01 mag and the
maximum magnitude difference is only 0.06 mag. In addition, the scatter of the
PL relation and the period–period ratio–luminosity relation is the same. Since
the metallicity is more dependent on period, we ultimately simplify the
period–period ratio–luminosity relation to the PL relation (see equation 3).
Here, we adopt the Wesenheit magnitude (Madore, 1982) i.e.
$W_{VI}=I-1.55(V-I)$ and $W_{G,BP,RP}=G-1.90(BP-RP)$.

$\displaystyle M_{W_{VI}}=(-4.523\pm 0.156)\log P_{\rm F}+(16.620\pm
0.048)-{\rm DM_{LMC}},\sigma=0.132\,{\rm mag},$ (3) $\displaystyle
M_{W_{VI}}=(-4.434\pm 0.153)\log P_{\rm 1O}+(16.079\pm 0.067)-{\rm
DM_{LMC}},\sigma=0.132\,{\rm mag},$ $\displaystyle M_{W_{G,BP,RP}}=(-3.623\pm
0.229)\log P_{\rm F}+(17.042\pm 0.071)-{\rm DM_{LMC}},\sigma=0.159\,{\rm
mag},$ $\displaystyle M_{W_{G,BP,RP}}=(-3.557\pm 0.225)\log P_{\rm
1O}+(16.606\pm 0.098)-{\rm DM_{LMC}},\sigma=0.159\,{\rm mag}$

We combine the LMC distance (Pietrzyński et al., 2019) and Gaia parallaxes to
calibrate the zero points of the PL relations as $M_{W_{G,BP,RP}}(P_{\rm
1O}=0.37{\rm d})=-0.348\pm 0.022$ mag and $M_{W_{VI}}(P_{\rm 1O}=0.37{\rm
d})=-0.496\pm 0.022$ mag. Meanwhile, we derive the zero point of the Gaia DR3
parallax, which is $13.4\pm 5.8\ \mu{\rm as}$. This parallax zero point agrees
well with that determined by classical Cepheids (Riess et al., 2021).

##  4 Distance and Metallicity Measurements

Figure 2: Metallicity and distance modulus (DM) determinations based on RRd
stars for M15 (red dots), IC 4499 (black dots), Sculptor (green dots), M3
(magenta dots), and M68 (blue dots). The mean values of RRd stars in five
targets are shown as open circles with $1\sigma$ internal errors. Solid
triangles indicate the parameters of the four globular clusters from the
Harris’s globular cluster catalog. To make the figure clearer, the distance
modulus of M15, IC 4499, and Sculptor are reduced by 0.5 mag, 0.8 mag, and 3.2
mag, respectively.

The established PL relation and PZ relation of RRd stars can be used to
determine the distances and metallicities of distant objects. We measure the
distances and metallicities of four globular clusters (IC 4499, M15, M3, and
M68) and a dwarf galaxy (Sculptor) based on RRd stars with Gaia data. Figure 2
shows the derived metallicities and distance moduli. The different colored
dots indicate individual metallicity and distance modulus determinations based
on RRd stars in M15 (red), IC 4499 (black), Sculptor (green), M3 (magenta),
and M68 (blue). Solid triangles indicate the results for the four globular
clusters from the globular cluster catalog of Harris (1996). The mean values
of RRd stars in five targets are shown as open circles with $1\sigma$ internal
errors. For dwarf galaxy Sculptor, the distance accuracy is about 1.8%, and
the metallicity dispersion is 0.2 dex. For globular clusters, the distance
accuracy is about 2-3%, and the metallicity dispersion is 0.05 dex. Recently,
we use RRd stars to measure the distances and metallicities of other dwarf
galaxies. For example, we measure the distances and metallicities of five
dwarf galaxies in M31, where the distance accuracy is about 1.8%. We also
determine the distance of the Fornax dwarf galaxy with an accuracy of 1.2%
based on $\sim 130$ RRd stars.

There are several advantages of adopting RRd stars as distance indicators.
Since the RRd star PL relation is not affected by metallicity, the galaxy
distance error can be optimized to 1.2% based on $\sim 100$ RRd stars. This
accuracy is as good as the distance accuracy of the LMC. RRd stars are not
rare, and they account for 3-5% of RR Lyr stars. The number of RRd stars in
the SMC, LMC, and the Galactic halo is comparable to the total number of
classical Cepheids in these galaxies. Moreover, distance measurements based on
RRd stars are less affected by extinction and binarity than the classical
Cepheids. RRd stars are the best distance indicators for dwarf galaxies. They
can also provide metallicity information for distance measurements of RRab
stars. In the future, the China Space Station Telescope and the Vera C. Rubin
Observatory Legacy Survey of Space and Time will discover more RRd stars. The
RRd star can be an independent distance ladder in the Local Group, used to
check the Cepheid (Riess et al., 2022) and the tip of the red giant branch
(Freedman et al., 2019) distance ladders.

##  5 Summary

We discover a period–(period ratio)–metallicity relation for RRd stars and
find it can predict metallicity as accurately as the low-resolution spectra.
We propose that the PL relation of RRd stars is not affected by the
metallicity and derive the PL relation of RRd stars. The zero point of the PL
relation is determined with an error of 0.022 mag (1.0% error in the
distance). We find that RRd stars can anchor galaxy distances with an accuracy
of 1-2% (dwarf galaxies) and 2-3% (globular clusters). In the Future, RRd
stars will be established as an independent distance ladder in the Local
Group.

## Acknowledgements

This work is supported by the National Natural Science Foundation of China
(NSFC) through the projects 12003046, 12173047, 11903045, 12133002, and
11973001. S.W. and X.C. acknowledge support from the Youth Innovation
Promotion Association of the CAS (grant No. 2023065 and 2022055).

## References

  * Braga et al. (2022) Braga, V. F., Fiorentino, G., Bono, G., et al. 2022, Monthly Notices of the Royal Astronomical Society, 517, 5368
  * Catelan & Smith (2015) Catelan, M., & Smith, H. A. 2015, Pulsating Stars. (Wiley-VCH)
  * Chen et al. (2018) Chen, X., Wang, S., Deng, L., et al. 2018, ApJS, 237, 28
  * Chen et al. (2020) Chen, X., Wang, S., Deng, L., et al. 2020, ApJS, 249, 18
  * Chen et al. (2023) Chen, X., Zhang, J., Wang, S., & Deng, L. 2023, Nat. Astron., arXiv:2306.10708
  * Clementini et al. (2022) Clementini, G., et al., et al. 2022, arXiv, 2206.06278
  * Cui et al. (2012) Cui, X.-Q., et al. 2012, Res. Astron. Astrophys., 12, 1197
  * Eisenstein et al. (2011) Eisenstein, D. J., et al. 2011, Astron. J., 142, 72
  * Freedman et al. (2019) Freedman, W. L., et al. 2019, ApJ, 882, 34
  * Harris (1996) Harris, W. 1996, ApJ, 112, 1487
  * Heinze et al. (2018) Heinze, A. N. Tonry J. L., Denneau L., et al. 2018, ApJ, 156, 241
  * Jayasinghe et al. (2018) Jayasinghe, T., et al. 2018, Monthly Notices of the Royal Astronomical Society, 477, 3145
  * Kovtyukh et al. (2016) Kovtyukh, V., et al. 2016, Monthly Notices of the Royal Astronomical Society, 460, 2077
  * Madore (1982) Madore, B. F. 1982, ApJ, 253, 575
  * Marconi et al. (2015) Marconi, M., et al. 2015, ApJ, 808, 50
  * Nemec & Moskalik (2021) Nemec, J. M., & Moskalik, P. 2021, MNRAS, 507, 781
  * Petersen (1973) Petersen J. O. 1973, A&A, 27, 89
  * Pietrzyński et al. (2019) Pietrzyński, G., et al. 2019, Nature, 567, 200
  * Popielski et al. (2000) Popielski, B. L., Dziembowski, W. A., Cassisi, S. 2000, Acta Astron., 50, 491
  * Riess et al. (2021) Riess, A. G., et al. 2021, ApJL, 908, L6
  * Riess et al. (2022) Riess, A. G., et al. 2022, ApJL, 934, L7
  * Sesar et al. (2017) Sesar, B. et al. 2017, ApJ, 153, 204
  * Soszyński et al. (2009) Soszyński, I. et al. 2009, Acta Astron., 59, 1
  * Soszynśki et al. (2016) Soszynśki, I., et al. 2016, Acta Astron., 66, 131
  * Soszynśki et al. (2019) Soszynśki, I., et al. 2019, Acta Astron., 69, 321


This paper deals with the problem of learning the probabilities of causation of subpopulations given finite population data. The tight bounds of three basic probabilities of causation, the probability of necessity and sufficiency (PNS), the probability of sufficiency (PS), and the probability of necessity (PN), were derived by Tian and Pearl. However, obtaining the bounds for each subpopulation requires experimental and observational distributions of each subpopulation, which is usually impractical to estimate given finite population data. We propose a machine learning model that helps to learn the bounds of the probabilities of causation for subpopulations given finite population data. We further show by a simulated study that the machine learning model is able to learn the bounds of PNS for $32768$ subpopulations with only knowing roughly $500$ of them from the finite population data.

§ INTRODUCTION

The probability of causation is a crucial concept that belongs to the third ladder of causality defined by Pearl <cit.> and plays a significant role in modern decision-making. The applications include the areas of marketing, political science, and health science. For example, using a linear combination of the probabilities of causation, Li and Pearl defined the benefit function of the unit selection problem, which is considered to be a revolution of the traditional A/B test heuristics <cit.>. Personalized decision-making, for another example, has been demonstrated that it should consider the probabilities of causation <cit.>. The label of a machine learning algorithm should also combine the probabilities of causation terms to capture the counterfactual behavior of the desired task <cit.>.

Using structural causal model (SCM) <cit.>, Pearl <cit.> first defined three basic probabilities of causation (i.e., PNS, PN, and PS). These probabilities of causation were then bounded tightly by Tian and Pearl <cit.> using Balke's programming <cit.>. The theoretical proof of those bounds was provided at <cit.>. After that, researchers started to use covariate information and the causal structure to narrow the bounds of the above probabilities of causation <cit.>. Li and Pearl <cit.> also presented extended studies of the probabilities of causation with nonbinary treatment and effect.

All listed works above are focused on providing bounds of the probabilities of causation on a specific population and require the experimental and observational distributions of the population. If the population can be divided into subpopulations by some observed characteristics, then the estimations of the probabilities of causation for each subpopulation require the experimental and observational distributions of each subpopulation. Consider the following example: an online music provider wants to increase its user subscription rate by sending gifts to new subscribers. Therefore, based on customer characteristics such as income, age, and usage, the company wants to identify customers who are likely to subscribe if and only if they received the gift. The gifts are value earphones; thus, the music provider prefers that these gifts be made only to those identified customers. Therefore, the music provider wants to know which kind of characteristics contains more desired customers. The quantities that the music provider wanted are then the PNS of subpopulations for each set of customer characteristics. If we want to apply Tian-Pearl's PNS bounds <cit.>, we need the experimental and observational distributions of each set of customer characteristics. However, the data available are usually finite many for the whole population; therefore, 1) it is impractical to estimate each distribution because the number of subpopulations is large; 2) some of the subpopulations are very rare (i.e., they appeared with a low probability) or even with no data associated; 3) we are not able to forecast the new coming subpopulations.

In this work, we are considering binary treatment and effect and will focus on PNS. It is not hard to extend to the other probabilities of caution. We also classify individual behavior into four response types, labeled complier, always-taker, never-taker, and defier <cit.>. Compliers are individuals who would respond positively if encouraged and negatively if not encouraged (i.e., PNS is the fraction of compliers). Always-takers are individuals who always respond positively whether or not they are encouraged. Never-takers are individuals who always respond negatively whether or not they are encouraged. Defiers are individuals who would respond negatively if encouraged and positively if not encouraged. We assume that the response type of an individual is determined by his characteristics (observed or unobserved characteristics); therefore, we propose a machine learning framework that is able to output the lower and upper bounds of PNS for each subpopulation given finite population data (i.e., learn the relations between the characteristics and the response types).

§ PRELIMINARIES

Here, the definitions of three basic probabilities of causation are reviewed <cit.>. We follow the language of counterfactuals in structural causal model in <cit.>.

We use $Y_x=y$ to denote the basic counterfactual sentence “Variable $Y$ would have the value $y$, had $X$ been $x$”. If not specified, we use $y_x$ to denote $Y_x=y$, $X$ stands for treatment, and $Y$ stands for effect. The experimental distributions in this paper are those in the form of the causal effects, $P(y_x)$, and the observational distributions in this paper are those joint probability function $P(x,y)$.

Three basic probabilities of causation are defined as follows:

Let $X$ and $Y$ be two binary variables in a causal model $M$, let $x$ and $y$ stand for the propositions $X=true$ and $Y=true$, respectively, and $x'$ and $y'$ for their complements. The probability of necessity is defined as the expression <cit.>

\begin{eqnarray}
\text{PN} & \delequal & P(Y_{x'}=false|X=true,Y=true)\nonumber \\
& \delequal & P(y'_{x'}|x,y) \nonumber
\label{pn}
\end{eqnarray}

\begin{eqnarray}
\text{PS} \delequal P(y_x|y',x') \nonumber
\label{ps}
\end{eqnarray}

\begin{eqnarray}
\text{PNS}\delequal P(y_x,y'_{x'}) \nonumber
\label{pns}
\end{eqnarray}

PNS stands for the probability that $y$ would respond to $x$ both ways, and therefore measures both the sufficiency and necessity of $x$ to produce $y$.

The tight bounds of PNS, PN and PS derived by Tian and Pearl <cit.> are then in the following forms:

\begin{eqnarray}
\max \left \{
\begin{array}{cc}
0, \\
P(y_x) - P(y_{x'}), \\
P(y) - P(y_{x'}), \\
P(y_x) - P(y)\\
\end{array}
\right \}
\le \text{PNS}
\le \min \left \{
\begin{array}{cc}
P(y_x), \\
P(y'_{x'}), \\
P(x,y) + P(x',y'), \\
P(y_x) - P(y_{x'}) + P(x, y') + P(x', y)
\end{array}
\right \}
\label{pnsub}
\end{eqnarray}

\begin{eqnarray*}
\max \left \{
\begin{array}{cc}
0, \\
\frac{P(y)-P(y_{x'})}{P(x,y)}
\end{array}
\right \}
\le \text{PN} \le
\min \left \{
\begin{array}{cc}
1, \\
\frac{P(y'_{x'})-P(x',y')}{P(x,y)}
\end{array}
\right \}
\label{pnub}
\end{eqnarray*}

\begin{eqnarray*}
\max \left \{
\begin{array}{cc}
0, \\
\frac{P(y')-P(y'_{x})}{P(x',y')}
\end{array}
\right \}
\le \text{PS} \le
\min \left \{
\begin{array}{cc}
1, \\
\frac{P(y_{x})-P(x,y)}{P(x',y')}
\end{array}
\right \}
\label{psub}
\end{eqnarray*}

To obtain bounds for a specific subpopulation, defined by a set $C$ of characteristics, the expressions above should be modified by conditioning each term on $C=c$. Therefore, if the experimental and observational distributions are available for every subpopulation, we are able to estimate the probabilities of causation of every subpopulation (here, we reviewed Tian-Pearl bounds for PNS, PS, and PN; Li and Pearl provided bounds for all types of probabilities of causation in <cit.>). However, in practice, some subpopulations have no adequate data (due to finite population data) to estimate their experimental and observational distributions. In this paper, we propose a machine learning model that takes the bounds of the probabilities of causation (i.e., those bounds that have adequate data to estimate) as the label and provides the estimations of the probabilities of causation of all rest subpopulations.

§ CAUSAL MODEL

In order to verify the accuracy of the learned bounds of PNS, we must first understand the data generating process to have the true PNS value and its bounds. The model we are using is shown in Figure <ref> (the coefficients in SCM are randomly generated as in <cit.>, see the appendix for the detail), where $X$ is a binary treatment, $Y$ is a binary effect, and $Z$ is a set of $20$ independent binary features (say $Z_1,...,Z_{20}$). The structural equations are as follow (for simplicity reason, we let $x=1,x'=0$, and $y=1,y'=0$):

[->,>=stealth',node distance=2cm,
thick,main node/.style=circle,fill,inner sep=1.5pt]
[main node] (1) [label=above:$Z$];
[main node] (3) [below left =1cm of 1,label=left:$X$];
[main node] (4) [below right =1cm of 1,label=right:$Y$] ;
[every node/.style=font=]
(1) edge node  (3)
(1) edge node  (4)
(3) edge node  (4);

The Causal Model, where $X$, $Y$ and $Z$ are binary treatment, binary effect, and $20$ independent binary features, respectively.

\begin{eqnarray*}
Z_i &=& U_{Z_i} \text{ for } i \in \{1,...,20\},\\
%X &=& 1 \text{ if } M_X +U_X > 3,\\
X&=&f_X(M_X,U_X) =
\left \{
\begin{array}{cc}
1& \text{ if } M_X +U_X > 0.5, \\
0& \text{otherwize},
\end{array}
\right \}\\
%Y &=& 0.867991236792X+M_Y+U_Y,\\
Y&=&f_Y(X,M_Y,U_Y) =
\left \{
\begin{array}{cc}
1& \text{ if } 0< CX+M_Y+U_Y < 1\\
1& \text{ if } 1 < CX+M_Y+U_Y < 2,\\
% 1& \text{ if } 1 < CX+M_Y+U_Y < 2, \\
0& \text{otherwize},
\end{array}
\right \}\\
&&\text{where, } U_{Z_i}, U_X, U_Y \text{ are binary exogenous variables with Bernoulli distributions,}\\ &&\text{C is a constant, and }M_X,M_Y\text{ are linear combinatations of }Z_{1},...,Z_{20}.
\end{eqnarray*}

The value of $C,M_X,M_Y$ and the distributions of $U_X,U_Y,U_{Z_1},...,U_{Z_{20}}$ for the model are provided in the appendix.

§ DATA GENERATING PROCESS

Based on the model defined in the last section, there are $20$ binary features. We made $15$ of them observable and $5$ of them unobservable. All exogenous variables are also made unobservable. We then have $2^{15}$ observed subpopulations.

§.§ Informer Data

The informer data must know the actual bounds of PNS of each subpopulation for comparison purposes. From the structural equation given in last section, the value of $X$ is determined by $U_X$ and $M_x$ (denoted by $f_X(M_X,U_X)$) and the value of $Y$ is determined by $X$, $M_Y$, and $U_Y$ (denoted by $f_Y(X,M_Y,U_Y)$). If all $20$ binary features are observable, then for a particular features $z=(z_1,...,z_{20})$, $M_X$ and $M_Y$ are fixed (denoted by $M_X(z)$ and $M_Y(z)$), then the PNS, experimental distribution, and observational distribution of this set of features are

\begin{eqnarray*}
PNS(z) &=& P(Y=0_{X=0}, Y=1_{X=1}|z)\\
&=& P(U_Y=0)*T_0 + P(U_Y=1)*T_1, \\
\text{where}, &T_0& =
\left \{
\begin{array}{cc}
1& \text{ if } Y(0,M_Y(z),0)=0 \text{ and } Y(1,M_Y(z),0)=1, \\
0& \text{otherwize},
\end{array}
\right\},\\
&T_1& =
\left \{
\begin{array}{cc}
1& \text{ if } Y(0,M_Y(z),1)=0 \text{ and } Y(1,M_Y(z),1)=1, \\
0& \text{otherwize}
\end{array}
\right \}.
\end{eqnarray*}

\begin{eqnarray*}
&=& P(U_Y=0)*Y(X,M_Y(z),0) + P(U_Y=1)*Y(X,M_Y(z),1).
\end{eqnarray*}

\begin{eqnarray*}
&=& P(U_X=0)*P(U_Y=0)*Y(X(M_X(z),0),M_Y(z),0)+\\
&&P(U_X=0)*P(U_Y=1)*Y(X(M_X(z),0),M_Y(z),1)) +\\ &&P(U_X=1)*P(U_Y=0)*Y(X(M_X(z),1),M_Y(z),0)) +\\ &&P(U_X=1)*P(U_Y=1)*Y(X(M_X(z),1),M_Y(z),1)).
\end{eqnarray*}

We assumed $15$ of the features are observable (say $Z_1,...,Z_{15}$), which means each subpopulation $c=(z_1,...,z_{15})$ consists $32$ sets of $20$ binary features (say $s_{0}=(z_1,...,z_{15},0,0,0,0,0), s_{1}=(z_1,...,z_{15},0,0,0,0,1), s_{2}=(z_1,...,z_{15},0,0,0,1,0), ...,s_{31}=(z_1,...,z_{15},1,1,1,1,1)$), then we have the PNS, experimental distribution, and observational distribution of all observed subpopulations are as follow:

\begin{eqnarray*}
PNS(c) &=& P(Y=0_{X=0}, Y=1_{X=1}|c)\\
&=& P(s_{0})/P(c)PNS(s_{0})+P(s_{1})/P(c)PNS(s_{1})+\\
&=& P(Z_{16}=0)P(Z_{17}=0)P(Z_{18}=0)P(Z_{19}=0)P(Z_{20}=0)PNS(s_{0})+\\
\end{eqnarray*}

\begin{eqnarray*}
&=& P(Z_{16}=0)P(Z_{17}=0)P(Z_{18}=0)P(Z_{19}=0)P(Z_{20}=0)P(Y=1|do(X),s_{0})+\\
&& P(Z_{16}=0)P(Z_{17}=0)P(Z_{18}=0)P(Z_{19}=0)P(Z_{20}=1)P(Y=1|do(X),s_{1})+\\
&& P(Z_{16}=0)P(Z_{17}=0)P(Z_{18}=0)P(Z_{19}=1)P(Z_{20}=0)P(Y=1|do(X),s_{2})+...+\\
&& P(Z_{16}=1)P(Z_{17}=1)P(Z_{18}=1)P(Z_{19}=1)P(Z_{20}=1)P(Y=1|do(X),s_{31}).
\end{eqnarray*}

\begin{eqnarray*}
&=& P(Z_{16}=0)P(Z_{17}=0)P(Z_{18}=0)P(Z_{19}=0)P(Z_{20}=0)P(Y=1|X,s_{0})+\\
&& P(Z_{16}=0)P(Z_{17}=0)P(Z_{18}=0)P(Z_{19}=0)P(Z_{20}=1)P(Y=1|X,s_{1})+\\
&& P(Z_{16}=0)P(Z_{17}=0)P(Z_{18}=0)P(Z_{19}=1)P(Z_{20}=0)P(Y=1|X,s_{2})+...+\\
&& P(Z_{16}=1)P(Z_{17}=1)P(Z_{18}=1)P(Z_{19}=1)P(Z_{20}=1)P(Y=1|X,s_{31}).
\end{eqnarray*}

The informer view of the bounds of $PNS(c)$ (i.e., true bounds) could be obtained using Equation <ref> and above observational and experimental distributions.

§.§ Experimental Sample

Here is how we collected $5000000$ experimental samples for training purposes. From the causal model in Section <ref>, an individual is determined by $U_X$, $U_Y$, and $U_{Z_1},...,U_{Z_{20}}$. Therefore, we first randomly generated $(U_X,U_Y,U_{Z_1},...,U_{Z_{20}})$ using their distributions; then we randomly generated $X$ using $Bernoulli(0.5)$; the value of $Y$ is then $f_Y(X,M_Y,U_Y)$. We then obtained a experimental sample, $(U_{Z_1},U_{Z_2},...,U_{Z_{15}},X,Y)$ (since $15$ features are observable and identical to $U_{Z_i}$).

§.§ Observational Sample

Here is how we collected the $5000000$ observational samples for training purposes. Similarly to experimental data, we first randomly generated $(U_X,U_Y,U_{Z_1},...,U_{Z_{20}})$ using their distributions; the value of $X$ is then $f_X(M_X,U_X)$; the value of $Y$ is then $f_Y(X,M_Y,U_Y)$. We then collect a observational sample, $(U_{Z_1},U_{Z_2},...,U_{Z_{15}},X,Y)$ (since $15$ features are observable and identical to $U_{Z_i}$).

§ MACHINE LEARNING MODEL

§.§ Features and Label

The training features are $15$ observed features. We obtain the training label as follows: if a given set of features appeared more than $1300$ times (note we use the number $1300$ to have a precise estimation of PNS suggested by Li, Mao, and Pearl <cit.>) in those $5000000$ experimental samples and appeared more than $1300$ times in those $5000000$ observational samples, Frequentist will be used to estimate the experimental and observational distributions. We will then apply Equation <ref> to obtain lower bound and upper bound for this set of features (lower bound and upper bound are the labels). We totally have $529$ set of features of both lower and upper bounds for training purposes. The $529$ sets of features are then split into $423$ for the training set and $106$ for the testing set.

§.§ Learning

We used a simple fully-connected neural network to predict the lower and upper bounds from the feature set. Specifically, we use four multilayer perceptron (MLP) layers. ReLU is used for activation in the leading three layers while Sigmoid function is used for the output layer. We set the embeddings dimension as 128 for all layers, and train the model for 600 iterations with learning rate 0.01. Our experiments are done at an AWS p3.2xlarge instance.

§ EXPERIMENTAL RESULTS

We randomly selected $200$ subpopulations (among $32768$) and compared their learned bounds of PNS with the true bounds of PNS computed from the informer data. The results are shown in Figure <ref>. The learned bounds of PNS for subpopulations are a good fit for the true PNS bounds. The average error of the learned lower bound among $32768$ subpopulations is $0.0775$, and the average error of the learned upper bound among $32768$ subpopulations is $0.1371$. They are both acceptable errors given we have only $423$ training size to learn $32768$ subpopulations.

Lower bound of the PNS for subpopulations. These $200$ subpopulations are randomly selected from $32768$ subpopulations.

Upper bound of the PNS for subpopulations. These $200$ subpopulations are randomly selected from $32768$ subpopulations.

The learned bounds of PNS for subpopulations compared to the true bounds of PNS.

§ DISCUSSION

We demonstrated that the probabilities of causation for subpopulations could be learned from finite population data. However, we must discuss some properties of our proposed method further.

First, we applied the machine learning model to learn the bounds of PNS in this paper. This method can be applied to any probability of causation. Machine learning provides the ability to learn the probabilities of causation of subpopulations such that there is insufficient data to estimate the experimental and observational distributions. However, the key is that the label of the machine learning model should be the bounds of the probabilities of causation rather than experimental and observational data. Therefore, the size of the training set is no longer $5000000$ observational and experimental samples; it is the observational and experimental distributions and the bounds obtained from the distributions (in the example we had in the last section, our training data set are the size of $529$).

Second, we applied the simplest machine learning model. However, the key is the framework we proposed. The data generating process is also available; this is the first publicly available data generating process that can test approaches for counterfactual learning. Researchers who are familiar with fancy machine learning models are welcome to apply other machine learning models to this dataset.

Third, this work can be extended to unit selection problems <cit.> because the unit selection problems are a linear combination of the probabilities of causation. One can also apply this machine learning framework to predict the experimental distributions <cit.> using observational data.

§ CONCLUSION

We demonstrated how to obtain bounds of probabilities of causation for subpopulations using finite population data. We proposed a machine learning framework to deliver reasonable estimations. Experiments showed that the probabilities of causation defined by SCM are learnable with proper labels. Data-generating processes are also available for future machine learning models.

§ ACKNOWLEDGEMENTS

This research was supported in parts by grants from the National Science
Foundation [#IIS-2106908 and #IIS-2231798], Office of Naval Research [#N00014-21-1-2351], and Toyota Research Institute of North America

§ APPENDIX

§.§ The Causal Model

The first model in <cit.> are used where the coefficients for $M_X,M_Y$ and $C$ were uniformly generated from $[-1,1]$, and the Bernoulli distribution parameters were uniformly generated from $[0,1]$. The detailed model is as follows:

\begin{eqnarray*}
Z_i &=& U_{Z_i} \text{ for } i \in \{1,...,20\},\\
%X &=& 1 \text{ if } M_X +U_X > 3,\\
X&=&f_X(M_X,U_X) =
\left \{
\begin{array}{cc}
1& \text{ if } M_X +U_X > 0.5, \\
0& \text{otherwize},
\end{array}
\right \}\\
%Y &=& 0.867991236792X+M_Y+U_Y,\\
Y&=&f_Y(X,M_Y,U_Y) =
\left \{
\begin{array}{cc}
1& \text{ if } 0< CX+M_Y+U_Y < 1 \text{ or } 1 < CX+M_Y+U_Y < 2,\\
% 1& \text{ if } 1 < CX+M_Y+U_Y < 2, \\
0& \text{otherwize},
\end{array}
\right \}\\
&&\text{where, } U_{Z_i}, U_X, U_Y \text{ are binary exogenous variables with Bernoulli distributions.}
s.t., \\
&&U_{Z_1} \sim \text{Bernoulli}(0.352913861526), U_{Z_2} \sim \text{Bernoulli}(0.460995855543),\\
&&U_{Z_3} \sim \text{Bernoulli}(0.331702473392), U_{Z_4} \sim \text{Bernoulli}(0.885505026779),\\
&&U_{Z_5} \sim \text{Bernoulli}(0.017026872706), U_{Z_6} \sim \text{Bernoulli}(0.380772701708),\\
&&U_{Z_7} \sim \text{Bernoulli}(0.028092602705), U_{Z_8} \sim \text{Bernoulli}(0.220819399962),\\
&&U_{Z_9} \sim \text{Bernoulli}(0.617742227477), U_{Z_{10}} \sim \text{Bernoulli}(0.981975046713),\\
&&U_{Z_{11}} \sim \text{Bernoulli}(0.142042291381), U_{Z_{12}} \sim \text{Bernoulli}(0.833602592350),\\
&&U_{Z_{13}} \sim \text{Bernoulli}(0.882938907115), U_{Z_{14}} \sim \text{Bernoulli}(0.542143191999),\\
&&U_{Z_{15}} \sim \text{Bernoulli}(0.085023436884), U_{Z_{16}} \sim \text{Bernoulli}(0.645357252864),\\
&&U_{Z_{17}} \sim \text{Bernoulli}(0.863787135134), U_{Z_{18}} \sim \text{Bernoulli}(0.460539711624),\\
&&U_{Z_{19}} \sim \text{Bernoulli}(0.314014079207), U_{Z_{20}} \sim \text{Bernoulli}(0.685879388218),\\
&&U_{X} \sim \text{Bernoulli}(0.601680857267), U_{Y} \sim \text{Bernoulli}(0.497668975278),\\
\end{eqnarray*}

\begin{eqnarray*}
&M_X& =
\begin{bmatrix}
\end{bmatrix}\times
\begin{bmatrix}
0.259223510143\\ -0.658140989167\\ -0.75025831768\\ 0.162906462426\\ 0.652023463285\\ -0.0892939586541\\ 0.421469107769\\ -0.443129684766\\ 0.802624388789\\ -0.225740978499\\ 0.716621631717\\ 0.0650682260309\\ -0.220690334026\\ 0.156355773665\\ -0.50693672491\\ -0.707060278115\\ 0.418812816935\\ -0.0822118703986\\ 0.769299853833\\ -0.511585391002
\end{bmatrix},
M_Y =
\begin{bmatrix}
\end{bmatrix}\times
\begin{bmatrix}
-0.792867111918\\ 0.759967136147\\ 0.55437722369\\ 0.503970540409\\ -0.527187144651\\ 0.378619988091\\ 0.269255196301\\ 0.671597043594\\ 0.396010142274\\ 0.325228576643\\ 0.657808327574\\ 0.801655023993\\ 0.0907679484097\\ -0.0713852594543\\ -0.0691046005285\\ -0.222582013343\\ -0.848408031595\\ -0.584285069026\\ -0.324874831799\\ 0.625621583197
\end{bmatrix}
\end{eqnarray*}