Coincident detection significance in multimessenger astronomy

Introduction

On August 17th 2017, the observation by LIGO-Virgo of GW170817, a binary neutron star coalescence (BNS) (Abbott et al., 2017b, a), and by Fermi and INTEGRAL of GRB 170817A, a short gamma-ray burst (GRB) (Goldstein et al., 2017; Savchenko et al., 2017), began an unprecedented multimessenger observing campaign (Abbott et al., 2017c). Detections and non-detections across the electromagnetic (EM) spectrum and by neutrino observatories have already produced new insights and will continue to do so for some time yet.

Many of these insights critically depend on the significance of the association between the independent observations. Often, such significance is established by estimating a p-value, the probability of such an event or a more extreme event occurring under the null hypothesis that the observations originate from unrelated distinct sources. Specific applications include, e.g. Abbott et al. (2017a); Coulter et al. (2017); Soares-Santos et al. (2017) for GW170817 and its counterparts, Baret et al. (2012); Aartsen et al. (2014); Keivani et al. (2015) for offline triggered search methods, and Urban (2016) for online rapid identification. A small p-value demonstrates the data is inconsistent with the null hypothesis. The p-value cannot, though, be interpreted as the probability of the null hypothesis itself (Gelman et al., 2013). On the other hand, a large p-value does not necessarily imply that the null hypothesis has to be accepted, only that it cannot be rejected (Gregory, 2005).

We introduce a different, generic model comparison method to determine whether two events in separate data sets are produced by a common source or by unrelated phenomena. This Bayesian measure of significance asks fundamentally different questions compared to the Frequentist p-value approach: it quantifies a degree of belief or confidence when comparing two hypotheses, given a particular non-repeatable observation, while the p-value determines the consistency of the null hypothesis with the data and the error rate of determining significance (which is important for initial identification). (See Finn (1998) for a related discussion in the context of detection itself). The method is a direct comparison of the probabilities of alternative models and does not require empirical estimates of a background distribution for the interpretation of its result (although this may be necessary if the assumptions about the background are not trusted). Moreover, the framework requires explicit statements of the necessary assumptions; in particular, prior distributions on the relevant parameters and conditions for which significance can be factorized for different common model parameters (discussed later in Sec. 2.3). This approach is distinct from that of Kelley et al. (2013) in which the EM data is used as prior information to understand improvements in sensitivity for triggered searches.

In Sec. 2 we introduce the method in a general context; Eq. (16) is our primary result and describes how to calculate the Bayes factor for a common-source origin of two signals seen in separate data streams. In Sec. 3 we focus on the application of the method to multimessenger astronomy, considering a calculation of spatial and temporal significance. As an example, we apply it to the gravitational-wave and gamma-ray events GW170817 and GRB 170817A, showing that it strongly supports the hypothesis that they originate from a common source.

Generic derivation

Given two detections $a$ and $b$ in different data sets $D_{a}$ and $D_{b}$, we would like to assess the hypothesis that they originate from a common source. In general, the two detections will be described by different physical signal models $\mathcal{H}^{\textrm{S}}_{a}$ and $\mathcal{H}^{\textrm{S}}_{b}$, respectively. Each signal model will imply a likelihood, a set of parameters and an associated prior for those parameters. To assess whether they originate from a common source, the models must share a common set of parameters $\theta\in\Theta$.

We’ll use notation where $\mathcal{H}(\theta)\equiv[\mathcal{H}\;\textrm{and}\;\theta]$ denotes a hypothesis $\mathcal{H}$ with a particular choice for the parameters $\theta$, while $\mathcal{H}$ by itself denotes a hypothesis with unknown parameters, i.e. “for any choice of parameters $\theta$”. We can formally write this as $\mathcal{H}\equiv[\mathcal{H}(\theta)\text{ for any }\theta]$.

Then, we define the common-source hypothesis:

We also define $\mathcal{H}^{\textrm{N}}_{a/b}$ as the noise hypotheses (by which we mean any non-signal) for each data set. Then we can define any alternative hypothesis for which the observed detections in $a$ and $b$ are unrelated:

where $X$, $Y$ $\in\{\textrm{N},\textrm{S}\}$. We write this in a general form, but note that the noise hypothesis will not have any common model parameters. In total, there are four possible realizations of $\mathcal{H}^{XY}$, which we consider in detail below. However, $\mathcal{H}^{\textrm{S}\textrm{S}}$ is of particular interest in this work, being two unrelated signals from distinct sources.

These hypotheses imply priors on $\theta$ which in general differ from those implied by $\mathcal{H}^{\textrm{S}}_{a/b}$ individually: if a common source can only be detected in some subset of $\theta$, then $\mathcal{H}^{\textrm{C}}$ can only have prior support restricted to this subset. If this is not true, we identify the special case

The probability of the common-source hypothesis is given by

In this work, we will calculate the odds between $\mathcal{H}^{\textrm{C}}$ and different choices of $\mathcal{H}^{XY}$

is the Bayes factor and $P(\mathcal{H}^{\textrm{C}})/P(\mathcal{H}^{XY})$ is the prior odds. In Sec. 2.2 we discuss the calculation of the Bayes factor in general. The prior odds will depend on the context, but in Sec. 3.1.3 we calculate the prior odds modeling $\mathcal{H}^{\textrm{C}}$ and $\mathcal{H}^{\textrm{S}\textrm{S}}$ as realizations of a Poisson point process.

2 Derivation of the Bayes factor

If both data sets contain a signal from the same event, then they are not independent: $P(D_{a},D_{b}|\mathcal{H}^{\textrm{C}})\neq P(D_{a}|\mathcal{H}^{\textrm{C}})P(D_{b}|\mathcal{H}^{\textrm{C}})$. Instead we must compute

where the domain of the integral in the second line is restricted to the prior support of $\mathcal{H}^{\textrm{C}}$, namely

The need for this restriction arises because assuming that $(\theta=\theta^{\prime})$ and $\mathcal{H}^{\textrm{C}}$ are both true would be a contradiction if $P(\theta^{\prime}|\mathcal{H}^{\textrm{C}})=0$, and so $P(D|\theta^{\prime},\mathcal{H}^{\textrm{C}})$ would be undefined. Rearranging the likelihood in the integrand

where in the second step we have used that the likelihoods conditional on $\theta$ can be separated for the two data sets, provided that $\theta$ is the set of all model parameters common between the two likelihoods. In the last step, we again used that $P(\theta|\mathcal{H}^{\textrm{C}})>0$ within the integration interval. A subtle point is that $P(\theta|D_{a/b},\mathcal{H}^{\textrm{C}})$ is the posterior distribution for the common model parameters (given either $D_{a/b}$) marginalized over all other model parameters and using the prior implied by $\mathcal{H}^{\textrm{C}}$.

quantifies the agreement between the posterior distributions of $\theta$ derived independently. In this integral, the prior has the effect of setting a scale against which the degree of overlap can be compared.

Eq. (10)-(11) demonstrate how probabilities from separate data sets combine when each provides independent inferences about a common model parameter.

Returning to the Bayes factor, by our definition of the alternative hypothesis

We now specify three particular cases of interest for the alternative hypothesis. First, consider $\mathcal{H}^{\textrm{N}\textrm{N}}$: both $a$ and $b$ are caused by noise. Then Eq. (13) specializes to

where $\mathcal{B}_{\textrm{C}/\textrm{N}}$, in analogy with Eq. (6), is the common-source against noise Bayes factor. In the special case of Eq. (3), it can be shown that $\mathcal{B}_{\textrm{C}/\textrm{N}}(D_{a/b})=\mathcal{B}_{\textrm{S}/\textrm{N}}(D_{a/b})$, i.e. the independent signal against noise Bayes factor for each data stream.

This agrees with our intuition: if both signals are strong compared to the background noise and there is a good overlap of their common model parameters (quantified by the integral), we believe they originate from a common event. This is a powerful result as one can compute the joint Bayes factor from the common-source against noise Bayes factor for each detection individually, and the posterior overlap integral of $\theta$. Eq. (14) has analogous applications to the Fisher combined probability test used in Aartsen et al. (2014).

Second, consider $\mathcal{H}^{\textrm{S}\textrm{N}}$: $a$ was due to a signal, but $b$ was due to noise. For this case, Eq. (13) gives

For us to believe that detection $b$ is a real signal and originates from the same source as $a$, we require that the product of the Bayes factor for common-source against noise in $b$ and the posterior overlap be large. The case $\mathcal{B}_{\textrm{C}/\textrm{N}\textrm{S}}$ is analogous and the same special cases apply as mentioned previously.

Finally, consider $\mathcal{H}^{\textrm{S}\textrm{S}}$, the distinct-source hypothesis: both $a$ and $b$ are of the same nature as in the common-source hypothesis $\mathcal{H}^{\textrm{C}}$, but they are physically distinct (i.e. they belong to unrelated sources with different parameters $\theta_{a}\not=\theta_{b}$). Then,

This equation and the posterior overlap integral of Eq. (11) are the main results of this paper. This provides a simple and intuitive way to assess whether two detections originate from the same event, based on the posterior overlap of their common model parameters.

In the special case of Eq. (3), the prefactor to the posterior overlap integral is unity, such that

On the other hand, when Eq. (3) does not apply, the prefactor plays an important role in quantifying how the restricted prior implied by $\mathcal{H}^{\textrm{C}}$ affects the Bayes factor.

A similar result to Eq. (17) was obtained independently by Haris et al. (2017) in the context of strongly lensed gravitational wave signals from binary black hole mergers.

3 Factorization of the posterior overlap integral

When calculating $\mathcal{I}_{\theta}$, it is often convenient to factorize the posterior overlap integral, e.g. $\mathcal{I}_{\theta}=\mathcal{I}_{\phi}\mathcal{I}_{\psi}$ where $\phi\subsetneq\theta$ and $\psi=\theta\setminus\phi$. This factorization can only be performed, however, if $P(\phi|\psi,D_{A/B},\mathcal{H}^{\textrm{S}})=P(\phi|D_{A/B},\mathcal{H}^{\textrm{S}})$. There are situations in which this is the case, for example if the joint posterior distribution is an uncorrelated multivariate normal distribution. But generally, this will not be the case and the posterior over the full common parameter space must be used. There are however cases where, under certain assumptions, the integral can be approximately factorized. We will explore one such setting in the next Section.

Application to multimessenger transient astronomy

We now focus on the application of the above formalism to multimessenger transient astronomy. To guide our intuition, we consider a transient gravitational wave (GW) candidate and a detection made by an EM instrument, although, we could just as well consider any pair of EM, GW or neutrino detectors. Assuming that detections are made in both the GW and EM detectors and are independently significant, we aim to calculate $\mathcal{O}_{\textrm{C}/\textrm{S}\textrm{S}}(D_{\textrm{GW}},D_{\textrm{EM}})$, the odds quantifying the probability of the common-source hypothesis to a distinct-source hypothesis.

The Bayes factor should be calculated from all common source parameters; typically, this will involve parameters such as a characteristic time of the event, source direction, luminosity distance, and source orientation (Troja et al., 2017; Margutti et al., 2017). Ideally, the posterior overlap integral should be computed over the complete joint distribution of parameters since it will not generally factorize (see Sec. 2.3).

However, to illustrate the utility of the method, we will calculate the result considering only the spatial and temporal common parameters (specifically, the source direction $\boldsymbol{\Omega}$ and coalescence time of the BNS system $t_{\textrm{c}}$) and make assumptions under which the posterior overlap integral may be factorized. We also consider both observatories to be all-sky, neglecting non-isotropic and non-stationary sensitivity. The Bayes factor can then be calculated from Eq. (17) since the special case of Eq. (3) applies.

To calculate the Bayes factor, Eq. (17), we first write down the posterior overlap integral over the conditional joint distribution of the spatial and temporal parameters

We will now show that this can be factorized into a spatial and temporal overlap under the following assumptions. First, that the prior itself factors, $P(\boldsymbol{\Omega},t_{\textrm{c}}|\mathcal{H}^{\textrm{S}})=P(\boldsymbol{\Omega}|\mathcal{H}^{\textrm{S}})P(t_{\textrm{c}}|\mathcal{H}^{\textrm{S}})$. Second, that $t_{\textrm{c}}$ inferred from the GW data is exactly determined, i.e.

where the “hat” indicates the observed value. Then, Eq. (18) can be factorized as $\mathcal{I}_{\boldsymbol{\Omega},t_{\textrm{c}}}=\mathcal{I}_{t_{\textrm{c}}}\mathcal{I}_{\boldsymbol{\Omega}}$, where

This factorization is exact under the two assumptions made. However, the coalescence time is typically known with a nonzero uncertainty. For this case, taking $\widehat{t_{\textrm{c}}}$ to be a point estimate (the mean for example), the factorization is approximate, but applicable provided that over the uncertainty in $t_{\textrm{c}}$, $P(t_{\textrm{c}}|D_{\textrm{EM}},\mathcal{H}^{\textrm{S}}_{\textrm{EM}})$, $P(\boldsymbol{\Omega}|t_{\textrm{c}},D_{\textrm{GW}},\mathcal{H}^{\textrm{S}}_{\textrm{GW}})$, and $P(\boldsymbol{\Omega}|t_{\textrm{c}},D_{\textrm{EM}},\mathcal{H}^{\textrm{S}}_{\textrm{EM}})$ do not vary substantially. In Sec. 3.1.1 and Sec. 3.1.2, we will provide approximations for Eq. (20) and (21) under some reasonable assumptions and illustrate some of the subtleties in their calculation.

We note that a similar result to Eq. (21) was previously derived in Urban (2016); in particular, Eq. (3.6) of that work is equivalent to Eq. (21) assuming an isotropic prior. Then the resulting joint likelihood ratio is defined using the alternative hypothesis that went into Eq. (15).

To evaluate Eq. (20), the temporal overlap, we first need to consider how to compute $P(t_{\textrm{c}}|D_{\textrm{EM}},\mathcal{H}^{\textrm{S}}_{\textrm{EM}})$, the coalescence time given the EM observations. Typically, EM observations do not directly infer $t_{\textrm{c}}$, but some other well defined time $t_{\textrm{EM}}$, e.g. the time of peak luminosity. We therefore need to specify a model that relates these two times. One simple model is that both signals travel at the speed of light, but there is a delay $\Delta t=t_{\textrm{EM}}-t_{\textrm{c}}$ between the coalescence time and the EM emission which will depend on the physics (see, e.g. Finn et al. (1999); Abadie et al. (2012) for GRB delay time predictions), but also on how $t_{\textrm{EM}}$ is defined. To fold these predictions into the analysis, we must specify $P(\Delta t|\mathcal{H}^{\textrm{S}})$, a prior distribution on the delay-time (at the Earth), given the model. Assuming $\Delta t$ and $t_{\textrm{EM}}$ are independent, the posterior can be transformed as

where $p_{t_{\textrm{EM}}}(t_{\textrm{EM}})\equiv P(t_{\textrm{EM}}|D_{\textrm{EM}},\mathcal{H}^{\textrm{S}}_{\textrm{EM}})$ denotes the posterior distribution of $t_{\textrm{EM}}$.

Having defined how to relate the time inferred by the EM data to the coalescence time with a suitable model, we now calculate Eq. (20) under some simple assumptions. Eq. (19) was the first of these assumptions and was already applied in factorizing the full posterior overlap integral. In addition, let

Next we need a prior for the delay in the GW-EM arrival time, which could be due to differences in emission time or propagation speed of GW and EM radiation. For simplicity we take a uniform distribution,

That is, the EM emission can arrive any time between a minimum and maximum value compared to the GW-inferred coalescence time; outside of that interval, we are certain the two events are not related. Inserting these definitions into Eq. (22), we obtain

from which the numerator of Eq. (20) can be calculated.

1.2 Spatial overlap

We now discuss calculating Eq. (21), the spatial posterior overlap integral. The EM counterparts to GW events are expected to originate from the same source direction and hence $\mathcal{I}_{\boldsymbol{\Omega}}$ can be directly computed from Eq. (21).

To illustrate the subtleties of $\mathcal{I}_{\boldsymbol{\Omega}}$ and provide some intuition, in Fig. 1, we show four examples varying the size of the uncertainty region and angular separation of the means of the EM and GW sky localizations. For all examples, a uniform all-sky prior is used. In Fig. 1(a), the means of both posteriors are aligned, but the uncertainty on both is large with respect to the all-sky prior; therefore, $\mathcal{I}_{\boldsymbol{\Omega}}$ is greater than one, but not large enough to be of note. For Fig. 1(b), $\mathcal{I}_{\boldsymbol{\Omega}}$ strongly indicates the two detections are from the same event: the means are aligned and the uncertainties are small with respect to the all-sky prior. In Fig. 1(c) and (d), the means of the distributions are not aligned. While in (c) this results in modest evidence in favor of a common event, the separation is sufficiently wide in (d) to strongly disfavor a common source.

To help guide our intuition, we can also calculate $\mathcal{I}_{\boldsymbol{\Omega}}$ for the simplified case where the posterior distributions on the sky are uniform distributions, i.e. constant inside the sets $\Pi_{\textrm{GW}}$ and $\Pi_{\textrm{EM}}$ and zero outside. Labeling $\Delta(\Pi)$ the area of set $\Pi$ in square radians, we obtain

The $4\pi$ prefactor comes from the all-sky prior and acts as a metric to compare the size of the overlap. For example, if $\Pi_{\textrm{GW}}$ is entirely contained within $\Pi_{\textrm{EM}}$, then $\mathcal{I}_{\boldsymbol{\Omega}}(D_{\textrm{GW}},D_{\textrm{EM}})=4\pi/\Delta(\Pi_{\textrm{EM}})$: the Bayes factor is determined entirely by the fraction of the sky covered by the uncertainty on the EM detections (or vice versa if the EM posterior is contained within the GW posterior).

1.3 The spatial and temporal odds

This prior odds clearly depends on the co-observing time $T$. Combining this with the spatial and temporal Bayes factor (Eq. (21) and (28)) then gives

which is not dependent on the co-observing time. One special case is when $R_{\textrm{GW}}\simeq R_{\textrm{GW,EM}}\ll R_{\textrm{EM}}$, i.e. if signals detectable in EM only are much more frequent than in GW, but we otherwise have little information on the rates of GW detections with or without EM counterparts. This may typically occur if our estimates of $R_{\textrm{GW}}$ and $R_{\textrm{GW,EM}}$ are based on $\mathcal{O}(1)$ detection. The odds are then proportional to $1/(R_{\textrm{EM}}[\Delta t])$, which reproduces the temporal Bayes factor Eq. (28) setting $T=1/R_{\textrm{EM}}$, i.e. the waiting time between EM detections (where the great majority have no GW counterpart).

As can be expected intuitively, the association becomes less significant if the $\Delta t$ prior is broader or the prior background rate of signals is higher, but increases with the prior expected rate of joint detections.

2 Application to GW170817 and GRB 170817A

We now apply the example calculated in Sec. 3.1 to GW170817 and GRB 170817A, the result of which can be compared with Abbott et al. (2017a). We note that, the calculation presented here could be improved by using the full joint distribution without making assumptions that allow the result to be factorized, and including other pertinent model parameters such as the luminosity distance (for Fermi-GBM, this may be as simple as estimating the range of conceivable values).

The sky localization for the BNS inspiral and short GRB can be seen in Fig. 1 of Abbott et al. (2017a). Using the published localization FITS files (Singer, 2017; Goldstein et al., 2017) and a uniform prior distribution on the whole sky, Eq. (21) yields $\mathcal{I}_{\boldsymbol{\Omega}}=32.4$. The spatial overlap alone provides moderate support for the common-event model, the main limitation being the uncertainty on the localization of GRB 170817A.

These numbers are consistent with the p-values estimated in (Abbott et al., 2017a): the time overlap dominates, the spatial part is small but supports the hypothesis, and the overall factor is highly significant (the total p-value was found to be $5\times 10^{-8}$ (Abbott et al., 2017a)).

3 Comparison with p-values

There are parallels that can be drawn between the odds calculated in Sec. 3.1 and the p-value approach of (Abbott et al., 2017a); namely, the spatial overlap Eq. (21) with the $\mathcal{S}$ statistic and the form of the temporal overlap (i.e. inversely proportional to the background rate).

However, the two methods are not equivalent and the numerical values themselves cannot be directly compared, as they answer different questions. The odds is exactly our relative degrees of belief for the common- vs. distinct-source hypotheses, given the assumptions made in the calculation; the p-value tests whether the data is consistent or not with the null (distinct-source) hypothesis, and is typically interpreted as the rate at which a rule for deciding the significance of a joint detection leads to false positives (Finn et al., 1999).

Finally, a more practical difference is that while Eq. (21) can be directly interpreted as a Bayes factor for the spatial overlap, interpreting the $\mathcal{S}$ statistic requires numerical calculation of the background by randomly rotating a set of observed short GRB sky localizations.

Conclusions

We introduce a Bayesian model comparison approach to estimating our confidence that two multimessenger observations are due to a common source as opposed to an accidental coincidence of distinct sources. The primary result of this work is Eq. (16), which generically allows the calculation of the Bayes factor (and, hence, the odds) from the joint posterior distributions of common model parameters inferred independently from two data sets. This approach forces us to recognize the conditions under which the contributions to the Bayes factor can be factorized.

We provide an example where the spatial and temporal overlap calculation can be approximately factorized for two independent observations with isotropic observatories and apply the result to GW170817 and GRB 170817A. We find decisive evidence in favor of their association, consistent with Abbott et al. (2017a).

References