Get statistical significance from the likelihood and so from the posterior (in anisotropy expansion)

Question

I try to understand, from a technical point of view, how are computed the statistical significance from a Bayesian study (I guess) in this abstract below from article "Evidence for anisotropy of cosmic acceleration" by Jacques Colin, Roya Mohayaee, Mohamed Rameez, Subir Sarkar:

Observations reveal a 'bulk flow' in the local Universe which is faster and extends to much larger scales than is expected around a typical observer in the standard ΛCDM cosmology. This is expected to result in a scale-dependent dipolar modulation of the acceleration of the expansion rate inferred from observations of objects within the bulk flow. From a maximum-likelihood analysis of the Joint Lightcurve Analysis (JLA) catalogue of Type Ia supernovae we find that the deceleration parameter, in addition to a small monopole, indeed has a much bigger dipole component aligned with the CMB dipole which falls exponentially with redshift $z$: $q_0=q_m+\vec{q}_d\cdot \hat{n}\exp(−z/S)$. The best fit to data yields $q_d=−8.03$ and $S=0.0262$ ($⇒d∼100 \text{Mpc}$), rejecting isotropy ($q_d=0$) with $3.9\sigma$ statistical significance, while $q_m=−0.157$ and consistent with no acceleration ($q_m=0$) at $1.4\sigma$. Thus the cosmic acceleration deduced from supernovae may be an artefact of our being non-Copernican observers, rather than evidence for a dominant component of 'dark energy' in the Universe.

Indeed, I have few notions like the relation :

$$\text{posterior}=\frac{\text{likelihood}\times\text{prior}}{\text{evidence}}$$

using likelihood

or more classically :

$$p(\theta|d)=\frac{p(d|\theta)p(\theta)}{p(d)}$$

with $\theta$ are the parameters to estimate and $d$ are the data.

I would like to understand how the statistical significance announced (the first one $= 3.9 \sigma$ ) is computed from the Bayesian relations above.

I think this is computed from the posterior but how to get this value :

they estimate from the likelihood at $d_d = -8.03$ and $S = 0.0262$ : how to compute this $3.9 \sigma$ ? Do they use the MLE (Maximum Likelihood Estimator) or MAP (Maximum Aposteriori Probability) methods ?

I hope you will understand my issue of understanding since I am interested into the necessity to introduce a cosmological constant or not into standard model.

Any explanations are welcome.

Regards

Answer 1

FYI, you may get more insight into the statistical significance from two followup papers to your cited paper regarding the heated debate over the evidence of dark energy.

Here is the full story line:

The paper by Jacques Colin et al Evidence for anisotropy of cosmic acceleration attempted to falsify dark energy:

the cosmic acceleration deduced from supernovae may be an artefact of our being non-Copernican observers, rather than evidence for a dominant component of `dark energy' in the Universe.

A followup paper by David Rubin at el Is the expansion of the universe accelerating? All signs still point to yes falsified the above falsifying paper:

we outline the errors C19 makes in their treatment of the data and inference on cosmological parameters. Reproducing the C19 analysis with our proposed fixes, we find that the dipole parameters have little effect on the inferred cosmological parameters. We thus affirm the conclusion of RH16: the evidence for acceleration is secure.

In a rebuttal posted at arxiv this Monday A response to Rubin & Heitlauf: "Is the expansion of the universe accelerating? All signs still point to yes", the author of the original paper responded that:

We emphasize that our procedure is justified and that their criticism serves only to highlight the rather "arbitrary corrections" that are made to the data in order to infer isotropic cosmic acceleration. This is a vivid illustration of the 'Cosmological Fitting Problem' faced by observers who live in an inhomogeneous universe but still use the maximally symmetric FLRW cosmolgy to interpret observations.

I am expecting a rebuttal to the rebuttal to the rebuttal soonish.

Answer 2

As I understand it, all that is being done is to compare the calculated (log) likelihoods that have been calculated for the various models.

The best fit model has a higher log likelihood than a model with $q_d = 0$. How you go from a likelihood ratio to a significance of $3.9 \sigma$ depends on what criteria you choose. There are various Bayesian Information Criteria that can be used to assign a significance to how much one model is preferred over another. Which of these has been used is impossible to tell without a proper reference or link to the paper!

As far as I know it is not possible to simply reject a model based only on its maximum likelihood; only to compare models using their likelihood ratios.

Edit: I looked at the relevant paper - Jacques Colin et al. (2019). Their 3.9 sig,ma result is from a change in the "Bayesian Information Criterion" (see above link). To quote them:

Since ΔBIC between the model with qd = 0 and the model with qm = 0 is 9.86, this constitutes strong evidence against a universe that is accelerating isotropically.

In summary, I am not sure where the $3.9\sigma$ comes from, and the authors don't say.

EDIT: Clarification has been received by Rameez (would have been nice to see this stated in the paper). The authors have used Wilks' theorem, which asserts that the difference in $2\times$ the log likelihood between two competing (but nested) models is distributed as chi-squared, with a number of degrees of freedom given by the number of additional parameters in the more complex model.

The best-fitting model has $-2\ln L = -208.28$ (see Table 1 in the paper), and this compares with $-2 \ln L = -189.52$ for an isotropic model with no dipolar $q_d$. The difference in $2\ln L$ is 18.76 and there are two additional degrees of freedom (called $q_d$ and $S$). Checking a chi-squared table with 2 dof gives a p-value of 0.0000844. I calculate this to correspond to 3.76 sigma by integrating the normal distribution. I am unclear (given Rameez's clarification) why this differs from the 3.9 sigma quoted in the paper. The same calculation procedure equates p-values of 0.16 with 1-sigma, 0.0227 with 2-sigma, 0.00135 with 3 sigma and 1/3.5 million with 5 sigma (as per the often quoted Higgs boson discovery).

Answer 3

The 3.9 sigma comes from just Wilk's theorem. \Delta 2 LLH between 1st and second rows of table 1, compared against a Chi Square of 2 dof.

Answer 4

My coauthor Rameez says:

"The actual delta 2 llh is 18.767 (finite precision), corresponding to a p value of 0.0000792 according to scipy.stats.distributions.chi2 . Then if I do scipy.optimize.fsolve to find the sigma at which the 1 - normal distribution cdf becomes 0.0000792, I get 3.88 sigma. I must have rounded this to 3.9.

Scipy itself uses interpolations to do some of these so there could be finite precision issues there, though I can't find any and root gives me similar answers within 4 significant digits."

Hope this clarifies the issue. Thanks for your interest (and Rob Jeffries for the pedagogical explanation above).