An empirical approach to model selection: weak lensing and intrinsic alignments

The impact of gravitational lensing along a particular line of sight is determined by two quantities, known as convergence and shear. The convergence $\kappa$ term of the weak lensing transformation describes how much a galaxy on a particular line of sight is distorted due to intervening large scale structure. It is defined as the weighted mass overdensity $\delta$, integrated along the line of sight to the distance of the source $\chi_{\mathrm{s}}$:

where $\boldsymbol{\theta}$ is the angular position at which the source is observed. The kernel $W(\chi)$, defined in Eq. (5), is sensitive to the relative distances of the source and the lens. It is this geometrical term that makes cosmic shear sensitive to the expansion history of the Universe.

In general terms, the impact of intrinsic alignments (IAs) can be thought of as adding a coherent additional component to each galaxy’s shape. That is, the observed ellipticity is $\gamma^{\rm obs}=\gamma^{\rm G}+\gamma^{\rm I}+\epsilon_{SN}$, or the sum of a shear due to cosmological lensing, an IA-induced shear, and a random shape noise component. Although the latter is typically dominant for any single galaxy, it cancels when the ellipticity is averaged across a large population of galaxies. At the level of angular correlation functions, one has:

Again assuming the Limber approximation, the angular power spectra can be written as

Both the TATT and NLA pipelines include free parameters for redshift error and residual shear bias. We adopt the same modelling as Krause et al. (2021), giving us one $\Delta z$ and one $m$ parameter per redshift bin, or a total of eight nuisance parameters. Note however, that these parameters are prior dominated for Y3 cosmic shear-only chains, and so add relatively little to the effective dimensionality. For details about the priors see Appendix A and Table 2. We also adopt the fiducial DES Y3 cosmic shear scale cuts (see Krause et al. 2021 for an explanation of how these were derived).

For a given set of input parameters, we generate a noiseless DES Y3-like cosmic shear data vector, $\boldsymbol{D}$, using the theory pipeline described in Section 2. We assume the fiducial Y3 redshift distributions, as presented in Myles, Alarcon et al. (2021). In all data vectors, the same input flat $\Lambda$CDM cosmology is used $(\Omega_{\mathrm{m}}=0.29,\;A_{\mathrm{s}}=2.38\times 10^{-9},\;\Omega_{\mathrm{b}}=0.052,\;h=0.75,\;n_{\mathrm{s}}=0.99,\;\Omega_{\nu}h^{2}=0.00053)$. This corresponds to $\sigma_{8}=0.79$, $S_{8}=0.77$, where $S_{8}\equiv\sigma_{8}\sqrt{\Omega_{\mathrm{m}}/0.3}$. We choose these to match the marginalised mean values from the DES Y1 $3\times 2$pt chain used to generate IA samples (see Section 3.2 below). Note, however, that the exact values are not expected to affect our results. We also fix all the redshift and shear calibration nuisance parameters to zero when generating data vectors.

When constructing simulated data vectors, it is important to remember that IA model parameters are not independent. The total GI+II intrinsic alignment component in a scenario with, e.g., $A_{1}=A_{2}=1$ is very different from one with $A_{1}=0.1,A_{2}=1$. As a consequence, it is possible for two sets of input IA parameters to give similar cosmological bias (when analysed with NLA), but quite different $\chi^{2}$ values. Specific combinations of TATT parameter values may enhance or cancel out cosmological parameter bias, and so it is useful to sample the 5D TATT parameter space rather than scaling up individual parameters to explore the potential for cosmological parameter bias due to model insufficiency. Therefore, instead of a single mock data vector, we generate a set of 21 data vectors, all with the same cosmology, but with different possible IA scenarios. The number of mock data vectors is an analysis choice. The more we generate, the better we cover the IA parameter space, but it also increases computational costs. We verified that 21 was a sufficient number of scenarios to have samples presenting low, medium and high bias in cosmological parameters, while still being reasonable in terms of computational cost (i.e. the chains to run).

To do this, we follow the recipe set out in Section 2 of DeRose et al. (2019). Starting with the posterior from the DES Y1 $3\times 2$pt TATT analysis (the purple contours in Figure 1; Samuroff et al. 2019), we evaluate the covariance matrix of the $N_{p}=5$ TATT parameters, and perform an eigenvalue decomposition. We then use Latin Hypercube sampling to generate $N_{\rm samp}$ samples, which are roughly evenly distributed in $N_{p}$ dimensional space. Finally, we use the eigenvalues/vectors to rotate and normalise those samples into the parameter space. The results are shown in Figure 1. The idea is that this provides a slightly broader coverage than could be obtained simply by drawing points from the joint posterior, while maintaining the correlations between parameters. In this way, we can cover a range of marginal cases, which are pessimistic, but still consistent with the data; this is useful, since for our purposes it is more important to span the range of plausible TATT model parameters than to preserve the statistics of the Y1 posterior exactly.

Since real measurements unavoidably include an (unknown) noise realisation, the calibration of the bias-metric relation is inherently a probabilistic problem (we will return to this point in Section 5; see the discussion there for details). For this reason, it is important that our simulations capture all sources of scatter in the data.

For each of our 21 IA scenarios, defined by a set of input TATT parameter values $\boldsymbol{\theta}_{\mathrm{IA},i}$, we have a set of noisy data vectors $\tilde{\boldsymbol{D}}_{i,j}=\boldsymbol{D}(\boldsymbol{\theta}_{\mathrm{IA},i})+\boldsymbol{N}_{j}$, where noise realisation $\boldsymbol{N}_{j}$ is drawn from the covariance matrix, and is assumed to be independent of $\boldsymbol{\theta}_{\mathrm{IA},i}$. We use the final DES Y3 covariance matrix, which is analytic and includes a Gaussian shape noise and cosmic variance contribution, as well as higher order non-Gaussian and super-sample terms (Friedrich et al., 2021). In total we generate 50 noise realisations, which we apply to each data vector. This gives us a collection of 21 noiseless data vectors, and $21\times 50=1050$ noisy ones.

For testing, however, it is convenient to arbitrarily choose a single fixed noise realisation, which we refer to as fiducial noise. Figure 2 shows an example DES Y3-like data vector, generated using the setup described above, with the fiducial noise realisation added. For this particular example, the input TATT parameters are the mean values from the Y1 posteriors in Figure 1. For reference, we also show the noiseless version (purple solid), as well as the separate (again, noiseless) IA contributions. Since the IA signal in the lowest bin ($1,1$) seems to dominate, one could reasonably ask whether we could simply focus on this part of the data vector for model selection. We ultimately decide against this for a few reasons. First, although the IA signal is strongest in bin $1,1$, we can see there is also non-negligible signal in the surrounding bins (e.g. $2,1$ and $3,1$). Indeed, during the DES Y3 analysis, a range of IA mitigation techniques were explored, including dropping the lowest auto bin correlations. For moderate TATT scenarios, it was found that this could not reliably eliminate cosmology biases, suggesting the IA contamination is not confined to the $1,1$ data vector. Additionally, the degree to which the low redshift II contribution dominates the IA signal is somewhat dependent on the input TATT parameters, and so the strength of this assumption varies in parameter space.

Now that we have a set of noisy data vectors, the next step is to fit them with all parameters free. We analyse each data vector twice, fitting to our full set of cosmological and nuisance parameters, but in one case using TATT, and in the other NLA. We then define bias as the distance between the peak of the marginalised posteriors in the $S_{8}-\Omega_{\mathrm{m}}$ plane. Figure 3 illustrates this for a particular simulated data vector. In brief, the algorithm works by evaluating the Euclidean distance between the peaks of the two posteriors, $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ in $S_{8}-\Omega_{\mathrm{m}}$ space. It then sequentially computes the confidence ellipses of $\mathcal{P}_{1}$ at different $\sigma$ levels, and finds a value of $N_{\sigma}$ that minimises the distance between the ellipse and the peak of $\mathcal{P}_{2}$. Note that this is the same recipe used in Krause et al. (2021). The choice of $S_{8}$ and $\Omega_{\mathrm{m}}$ as the parameters of interest comes simply from the fact that these are the cosmological parameters best constrained by DES. One could conceivably use a more complicated separation metric that is sensitive to the full parameter space, along the lines of those used for assessing tensions between data sets (e.g. Lemos, Raveri, Campos et al. 2021). For our purposes, however, we follow Krause et al. (2021) and consider the simpler 2D metric to be sufficient.

We use noisy simulated data, as described in Section 3.3 – this is an important feature of our analysis, and it is necessary to allow us to meaningfully interpret our statistical metrics. Therefore, the relative separation of the two posteriors is a more useful quantity than the distance from the input values of $S_{8}$ and $\Omega_{\mathrm{m}}$. The $0.22\sigma$ value shown in Figure 3 is assessed relative to the TATT posterior. This found to be more stable than assessing it relative to the NLA posterior, particularly in relatively extreme IA scenarios where the NLA posteriors are significantly shifted and can be distorted by prior edge effects.

Finally, it is implicit in the above that marginalised TATT constraints represent correct results by which to measure bias. That is, when we refer to bias, we are in fact talking about bias in the cosmological model when assuming the NLA model, with respect to what we find when assuming the TATT model. Although this is clearly reasonable (since our data were created using TATT), marginal contours can be subject to projection effects. Indeed, since some of the TATT parameters are relatively poorly constrained by shear-shear analysis alone, the two IA models cannot be assumed to experience projection effects to the same degree. We test this in Appendix A, and find projection offsets between TATT and NLA at the level of $0.1\sigma$. This is well below the threshold of $0.3\sigma$ used for this work, and is thus unlikely to significantly affect our results.

We investigate two commonly used test statistics, the $\chi^{2}$ difference and the Bayes ratio. We show in Section 5 that the former is more robust against noise than the latter, and is therefore a more useful metric for the method we are proposing.

It is also worth briefly remarking that in our analysis we found about 50 (out of 1050) data points for which $\chi^{2}_{\rm TATT}>\chi^{2}_{\rm NLA}$, and so $\Delta\chi^{2}_{(\textrm{df})}<0$. Given that NLA and TATT are nested models, these points are unphysical (a more flexible model should always be able to produce a better or as-good fit). We conclude that they are an artefact of the sampling method; although we tested the robustness of our IS setup, and found it can reproduce the best fit from chains to reasonable precision, some level of sampling noise is still present. Given this logic, it is reasonable to assume that if we were to find $\Delta\chi^{2}_{\rm NLA}<\chi^{2}_{\rm TATT}$ in real data, there would likely be some follow-up investigation and the chains would be rerun. This is particularly true if the $\chi^{2}$ is an integral part of the analysis plan, as it is in our method. We thus choose to discard these points. It is worth bearing in mind, however, that this may or may not be a reasonable decision in other setups, depending on the models being compared and the details of the analysis.

Also note that, although our results are based on 21 IA samples, we did initially draw 25 scenarios (see Section 3.2). Of these 25, we found four to be so extreme that the NLA PolyChord chains failed to converge in the noiseless case. These resulted in highly distorted and often bimodal contours in the $S_{8}-\Omega_{\mathrm{m}}$ plane, making it difficult to obtain meaningful estimates for the bias. Given this, and also for the reasoning discussed above, we choose to omit these samples from our analysis. This leaves us with a total of 21 IA samples.

Given all the definitions set out in the sections above, we now follow the recipe outlined below, in order to map and calibrate the bias-metric relation. These steps are, in essence, our method; when written out in this form, it can be very easily generalised to other model selection problems beyond our particular example of IA in cosmic shear.

1. 1.

   Sample IA scenarios: Draw about $10-30$ parameter samples from either a posterior from a previous analysis or from some reasonable priors using the method described in Section 3.2 (we used 21 drawn from DES Y1 TATT posteriors).

2. 2.

   Generate data vectors: Generate a simulated noiseless data vector for each IA scenario drawn in the previous step. Other model parameters (e.g., cosmological and nuisance parameters) should be fixed to some fiducial values. See Section 3.1.

3. 3.

   Analyse noiseless data vectors: There are two chains per data vector, one corresponding to model $M_{1}$, and the other to $M_{2}$. Again, we compared the TATT and NLA IA models using PolyChord to compute statistics. These choices might vary under different applications. Details on the sampling can be found in Section 3.4.1.

4. 4.

   Compute parameter bias and plot out the bias-metric relation: Demonstrate that the noiseless data vectors show a clear correlation between the test statistic (e.g., $\Delta\chi^{2}_{(\textrm{df})}$ or Bayes factor) and parameter bias, as described in Section 4.1.

5. 5.

   Generate noise realisations: For each data vector, generate $\sim 50$ noise realisations using the covariance matrix, as explained in Section 3.3. In our case, that gives a total of $21\times 50=1050$ noisy data vectors.

6. 6.

   Analyse noisy data vectors: As discussed above, we choose to use importance sampling to give a fast approximation for the noisy posteriors. For all noise realisations (50) and IA scenarios (21), for both IA models (2) $\longrightarrow$ 2100 total, estimate the posterior, the NLA-TATT bias and the model test statistics. See Section 3.4.2.

7. 7.

   Calculate probability: Plot out the bias-metric relation. Use the quantities computed in the previous step to calculate the probability of bias greater than some pre-defined threshold $X\sigma$ (see Section 4.2.3).

8. 8.

   Run analyses on blinded data: Run a full chain on the real data in order to evaluate the observed model statistic. This can then be interpreted in terms of bias probability using the results of the above step. See Section 5.

Note that, to obtain an accurate calibration of the $\Delta\chi^{2}_{(\textrm{df})}$ value, all aspects of the modelling should be as close to the final fiducial analysis setup as possible. For an estimate of the computational resources required to employ the proposed method, see Appendix E.

Considering first the noiseless case, Figure 4 shows the relation between bias in the $S_{8}-\Omega_{\mathrm{m}}$ plane and the NLA-TATT $\Delta\chi^{2}_{(\textrm{df})}$. Each point results from running two chains on the same noiseless simulated data vector, first using NLA, and then using TATT. As defined in Eq. (20), large values of $\Delta\chi^{2}_{(\textrm{df})}$ indicate statistical preference for the larger model (i.e., TATT). We see a relatively tight relation between bias and $\Delta\chi^{2}_{(\textrm{df})}$, going from $\Delta\chi^{2}_{(\textrm{df})}<1$ when bias is small to relatively large values at the high bias end: low bias $\rightarrow$ small $\chi^{2}$ difference, high bias $\rightarrow$ large $\chi^{2}$ difference. Interestingly, the relation appears to have the form (approximately) of a double power law, with a steep slope in the high bias regime, switching to a somewhat shallower function below $0.3\sigma$. It is worth stressing, however, that this relation is empirical. We do not have a particular expectation for its shape, and it is likely that the details depend on the analysis choices and survey setup. Note that even without data vector noise, this relation presents some scatter. This arises both from sampler noise, and from the fact that this is a complex high dimensional problem, for which two sets of IA values that produce biases of a similar magnitude will not necessarily produce identical $\Delta\chi^{2}_{(\textrm{df})}$ values. The vertical dotted line marks the $0.3\sigma$ bias threshold adopted by DES (Krause et al., 2021), which we adopt as our fiducial tolerance (see Section 5.2 below). Although we cannot use this noiseless result for any empirical method because real data will always contain noise, confirming that these quantities clearly correlate is a necessary first step in our method, and important to check before incurring the expense of further calculations. We will see in the next section that the correlation between $\Delta\chi^{2}$ and $S_{8}-\Omega_{\mathrm{m}}$ bias holds (with some additional scatter) when we proceed to the noisy case.

In addition to the $\Delta\chi^{2}_{(\textrm{df})}$, we also consider the Bayes Ratio as a potential model comparison metric; while the former presents a clear relation to the bias (as seen in Figure 4), we find the latter be a relatively weak indicator, with additional intrinsic scatter. This can be seen in Appendix C, and in particular Figure 11. We also note that $R$ and $\Delta\chi^{2}_{(\textrm{df})}$ are correlated. In principle one could seek to combine them, but naively treating them as independent metrics is almost certainly double-counting information. Therefore, here we focus on the results using $\Delta\chi^{2}_{(\textrm{df})}$. For further discussion of the Bayes ratio see Appendix C.

In Figure 5 we illustrate the impact of data vector noise in the bias-$\Delta\chi^{2}_{(\textrm{df})}$ plane. We show the same 21 noiseless samples discussed above, but now overlain with multiple different noise realisations, as approximated using importance sampling. As we can see, noise introduces scatter in the bias-$\Delta\chi^{2}_{(\textrm{df})}$ relation. While this noise is considerably less than in the case of the Bayes factor (for which we show in Appendix C that the scatter due to noise is so large that the relation with bias is extremely weak), it is still non-negligible.

For comparison, we show the $\Delta\chi^{2}_{(\textrm{df})}$ cut-offs implied by some standard model selection metrics: BIC, AIC and a $p-$value significance threshold⁴⁴4Assuming that the TATT model has 2 additional effective degrees of freedom compared with NLA. Note that this $df_{\text{diff}}$ value was calculated for the DES Y3 shear-only (no shear ratio) case in Secco, Samuroff et al. (2022). It is thus valid for our particular case, but would not necessarily hold under changes to the data vector or analysis choices. of $p=0.05$ (see Section 4.2.1 for definitions). Unfortunately, in the presence of noise, we see that all three cut-offs are relatively weak indicators of bias – i.e., they still favour the simpler model even when significant amount of bias has been introduced in the cosmological parameters. Even in the case of AIC, which is the strictest of the three, there are a problematic fraction of noise realisations where the observable metric favours NLA, and yet NLA is biased by $>0.3\sigma$ (see the points in the lower right hand corner of Figure 5). This illustrates a key motivation for adopting an empirical calibration. Theoretical limits imposed using, e.g., $p-$values are not designed to optimise the quantities we care most about (i.e., parameter biases). For a given analysis, it is impossible to know from first principles what level of bias is excluded for a given statistical metric cut-off without some form of calibration.

These observations have important consequences for our method. Since the exact noise realisation in any real data set is unknown, one cannot simply run a single set of IA samples (as in Figure 4), and perform a 1:1 bias-$\Delta\chi^{2}_{(\textrm{df})}$ mapping. Nor, as we can see from Figure 5, can we simply fall back on theoretical cut-offs to reliably guard against model bias. Instead, we must consider the problem as a probabilistic one, and factor in the uncertainty from noise.

It is also worth taking a moment to consider a related question: if the true IA scenario is extreme enough to give significant cosmological biases, would there be clear red flags from NLA alone, assuming that no fits were carried out with TATT? If this were the case, it would provide a simpler route – instead of performing empirical calibration using simulations, we could simply run one model on the data, and interpret results to see if a more sophisticated model is needed. Considering our 21 IA scenarios with a fixed noise realisation, we find that around $50\%$ of cases with bias $>0.3\sigma$ have $\chi^{2}$ values that appear entirely consistent with being drawn from the corresponding $\chi^{2}$ distribution⁵⁵5Where the NLA $p-$value here is calculated by assuming $\chi^{2}_{\rm NLA}$ is drawn from a $\chi^{2}$ distribution with 224 degrees of freedom (see Secco, Samuroff et al. 2022). The null hypothesis in this case is that the NLA model is adequate to describe the data, and so small $p-$values would indicate model insufficiency., $p(\chi^{2}_{\rm NLA})>0.05$. A similar picture is seen when we consider a single high bias IA scenario with alternative noise realisations – computing $p-$values for each realisation, the majority are above 0.05, even in the presence of bias $>0.3\sigma$. In other words, even in cases where NLA is significantly biased, it is not necessarily obvious from considering the uncalibrated value of $\chi^{2}_{\rm NLA}$ alone. In contrast, the method we propose, using $\Delta\chi^{2}_{(\textrm{df})}$, can correctly identify the need to use TATT to achieve sufficiently unbiased results $77.1/77.5\sim 99\%$ of the time (see the confusion matrix in Section 5.2).

Likewise, although extreme biases do tend to distort the shape of the posteriors, this is not always true in more moderate (but still significantly biased) cases. Figure 8 shows the NLA posteriors in a few different IA scenarios, spanning the range from almost no bias (purple shaded), to $\sim 1\sigma$ (pink, open contours). For reference we also show the posteriors from TATT fits to the same data vectors in Appendix B. Taken in isolation, none of these show clear signs of problems with the model. It is also interesting that IA mismodelling bias does not always translate into significantly non-zero values for the inferred NLA parameters. In the medium bias case, for example, $A_{1}$ and $\eta_{1}$ are both consistent with zero to $<1\sigma$. Here there is a relatively strong degeneracy between $A_{1}$ and $S_{8}$, allowing both $A_{1}\sim 0$ combined with low $S_{8}$, but also a stronger IA amplitude $(A_{1}\sim 1)$ with a larger $S_{8}$. In projection, this results in broad contours on both parameters (notice the black contours in the upper panel of Figure 8 are slightly wider than the others, with a more prominent asymmetry at high/low $S_{8}$).

The results discussed so far have a number of direct implications for the question of intrinsic alignment model selection. Primarily, we have shown that it is possible to perform empirical model selection with lensing data. There is a clear relation between cosmological parameter bias and $\Delta\chi^{2}_{(\textrm{df})}$, which allows one to define a threshold that can then be applied to the real data. That said, the failure of conventional statistical metrics (e.g., $p-$values) to identify scenarios with significant cosmological parameter biases is notable, and should be kept in mind when trying to understand statistics derived from any single run on real data. The properly calibrated model statistics, however, provide an alternative to the model selection exercises used in previous analyses, which have tended to rely on either simulated analyses (Secco, Samuroff et al. 2022 Section A3), or arguments based on direct-detection studies (Hikage et al. 2019; Joachimi et al. 2021a, Sections 5.4 and 2.4 respectively). The empirical method is arguably an advance on both; first of all, it avoids questions about what constitutes an “extreme" model, which tend to arise in the simulation-based approach. Since the current best constraints on TATT model parameters are relatively weak, it is relatively easy to select an IA scenario that is both consistent with observations, and which would cause significant bias in an NLA analysis (note that this is still true in light of the most recent DES Y3 results Secco, Samuroff et al. 2022; Amon et al. 2022; DES Collaboration 2022). Our empirical approach also avoids the uncertainties that are inherent in extrapolating observations based on direct IA measurements on one very specific type of galaxies to weak lensing measurements on another population entirely.

Although the empirical approach has various strengths, is worth reiterating that data vector noise is a significant source of scatter in the bias-$\Delta\chi^{2}_{(\textrm{df})}$ relation. For this reason it is important to accurately simulate the noise properties of the particular data set. While this is in principle simple, given an accurate covariance matrix, it does mean the model selection exercise needs to be repeated for any new data set or changes to the analysis. It also means it is crucial to have a fast and accurate way of estimating posteriors for a large number of noise realisations, such as the IS framework used here.

It is also interesting, finally, to consider how our findings relate to the real Y3 results. Although comparing with the full $3\times 2$pt results is difficult, for the reasons given above, Secco, Samuroff et al. (2022) (Section VIIB and Table III) present a comparison of IA models without shear ratios, an analysis configuration that matches ours. Specifically, comparing the 2 parameter NLA model with 5 parameter TATT, they find $R=1.70\pm 0.36$ and $\Delta\chi^{2}_{(\textrm{df})}=5.2/2=2.6$. Interpreted with the help of Figure 6, this puts the risk of NLA being biased by $>0.3\sigma$ at somewhere around $30\%$, meaning runs using NLA on Y3 were more likely than not to be unbiased to within the $0.3\sigma$ threshold.