3.1. Basics

Models for the data are built on assumptions about the nature of both the noise and signal being measured. The data model is described by the likelihood function

(1) \begin{align} {\mathcal{L}}(d | \theta) ,\end{align}

where d is the data and $\theta$ is a set of parameters describing the noise and/or signal.Footnote b The likelihood function is a normalised probability density function for the data, not for the parameters $\theta$ (see Thrane & Talbot Reference Thrane and Talbot2019, for more details):

(2) \begin{align} \int d(d) \, {\mathcal{L}}(d| \theta) = & 1 , \end{align}

(3) \begin{align} \int d\theta \, {\mathcal{L}}(d| \theta) \neq & 1 .\end{align}

It is useful to define a marginal likelihood, which is also known as the Bayesian evidence:

(4) \begin{align} {\mathcal{L}}(d) = \int d\theta {\mathcal{L}}(d|\theta) \pi(\theta) .\end{align}

Here, $\pi(\theta)$ is the prior distribution for the parameters $\theta$ . The marginal likelihood is a model for the data, averaged over realisations of $\theta$ .

A well-known example of a model for gravitational-wave data is the Whittle likelihood model for Gaussian time-series noise, shown here for a single frequency bin:

(5) \begin{align} {\mathcal{L}}(\tilde{d}|\theta) = \frac{1}{2\pi\sigma^2} e^{-|\tilde{d}|^2 / 2\sigma^2} .\end{align}

Here, $\tilde{d}$ represents the frequency-domain gravitational-wave strain while $\sigma^2$ is related to the noise power spectral density (PSD) P and the frequency bin width $\Delta f$ :

(6) \begin{align} \sigma^2 = \frac{P}{4\Delta f} .\end{align}

Equation (5) is an example of a parameter-free model of the data. If, for example, a gravitational-wave signal from a compact binary coalescence is present, then the likelihood depends on the ${\gtrsim}15$ parameters associated with a compact binary coalescence (component masses, spins, etc.),

(7) \begin{align} {\mathcal{L}}(\tilde{d}|\theta) = \prod_k \frac{1}{2\pi\sigma_k^2} e^{-\left|\tilde{d}_k-\tilde{h}_k(\theta)\right|^2 / 2\sigma_k^2} .\end{align}

Here, $\tilde{h}_k(\theta)$ is a model for the frequency-domain strain from a gravitational-wave signal given binary parameters $\theta$ in frequency bin k. Sometimes, $h_k(\theta)$ , which can be defined in either the time domain or the frequency domain, is referred to as ‘the waveform model.’

Examining Equation (7), we can see various ways in which the likelihood can be misspecified, which we represent diagrammatically in Figure 2. First, the waveform $h(\theta)$ may be misspecified, which can lead to well-documented systematic error (Ohme et al. Reference Ohme, Nielsen, Keppel and Lundgren2013; Wade et al. Reference Wade, Creighton, Ochsner, Lackey, Farr, Littenberg and Raymond2014; Ashton & Khan Reference Ashton and Khan2020; Gamba et al. Reference Gamba, Breschi, Bernuzzi, Agathos and Nagar2021; Huang et al. Reference Huang, Haster, Vitale, Varma, Foucart and Biscoveanu2021). This is an example of a misspecified signal model. Second, the noise model can be misspecified. There are typically two ways that this can happen. One possibility is that the functional form of the likelihood is correct, but the noise PSD is misspecified. A number of papers have proposed different methods for estimating the noise PSD in order to minimise this form of misspecification, for example, Littenberg & Cornish (Reference Littenberg and Cornish2015), Cornish & Littenberg (Reference Cornish and Littenberg2015), Chatziioannou et al. (Reference Chatziioannou, Haster, Littenberg, Farr, Ghonge, Millhouse, Clark and Cornish2019).

Figure 2. Forms of misspecification that we explore in this Article. Individual events can be misspecified if the model for the noise or the signal is not an adequate description of reality. The population of events may also be misspecified. This manifests itself as prior misspecification, which can impact both individual analyses (where the prior may be restricted to a limited portion of the true extent of the posterior) and population analyses (where the goal is to uncover the true distribution of the population).

The other possibility is that the functional form of the likelihood is itself misspecified. This may occur because of non-Gaussian noise artefacts (Röver, Meyer, & Christensen Reference Röver, Meyer and Christensen2010) or uncertainty in the noise PSD (Talbot & Thrane Reference Talbot and Thrane2020; Biscoveanu et al. Reference Biscoveanu, Haster, Vitale and Davies2020; Banagiri et al. Reference Banagiri, Coughlin, Clark, Lasky, Bizouard, Talbot, Thrane and Mandic2020), both of which yield broader tails than the Whittle distribution. Likewise, marginalising over calibration uncertainty broadens the likelihood function (Sun et al. Reference Sun2020; Payne et al. Reference Payne, Talbot, Lasky, Thrane and Kissel2020; Vitale et al. Reference Vitale, Haster, Sun, Farr, Goetz, Kissel and Cahillane2021).Footnote c Covariance between neighbouring frequency bins induced from finite measurements of continuous noise processes can also lead to misspecification if not correctly accounted for (Talbot et al. Reference Talbot, Thrane, Biscoveanu and Smith2021; Isi & Farr Reference Isi and Farr2021). In the subsequent Subsections, we describe tests for these different forms of misspecification.

3.2. Testing for a misspecified signal model

In order to test for a misspecified waveform, it is sometimes useful to look at the whitenedFootnote d residuals of the data in frequency

(8) \begin{align} \tilde{r}(f|\theta) \equiv \frac{\tilde{d}(f)-\tilde{h}(f|\theta)}{\sigma(f)} ,\end{align}

and time

(9) \begin{align} r(t|\theta) = {\mathcal{F}}^{-1}\left[\tilde{r}(f|\theta)\right] ,\end{align}

where ${\mathcal{F}}^{-1}$ is the discrete inverse Fourier transform. Residuals can be useful to test for waveform misspecification because the differences between waveform models are clearly seen in the time and frequency domain.Footnote e Additionally, if there is a terrestrial noise artefact (glitch) present in the data, it is likely to appear clearly in the residuals. For example, in Abbott et al. (Reference Abbott2016b), the best-fit, time-domain residuals for GW150914 were shown to be consistent with Gaussian noise (see Figure 1 of Abbott et al. Reference Abbott2016b), helping to show that the data are well explained by a gravitational waveform in Gaussian noise.

Recommended tests:

• • Plot the whitened residuals in both time r(t) and frequency $|\tilde{r}(f)|$ and visually inspect for consistency with zero. In the frequency domain, the set of $\tilde{r}(f)$ are approximately independent measurements with Gaussian uncertainty equal to unity. Include residual curves for many posterior-distribution draws of $\theta$ in order to show theoretical uncertainty. A word of caution: the time-domain residuals are highly covariant in real data, and so it is less straightforward to interpret misspecification in the time domain than in the frequency domain.

• • As a first step, look for consistency by eye. If there are signs of misspecification, one may quantify the inconsistency, for example, ‘the residuals are five standard deviations away from zero at $500\ \mathrm{Hz}$ .’ While post hoc analysis is helpful for catching egregious misspecification, in an ideal world, one should define the consistency tests a priori for unbiased tests. In practice, this is not always possible.

Demonstration: Footnote f Our signal model is a sine-Gaussian chirplet (i.e., a sine wave multiplied by a Gaussian function). We create two synthetic sets of data with Gaussian noise. The correctly specified data contains a signal that matches our model. The second dataset contains an intentionally misspecified signal: the same sine wave as before, but multiplied by a Tukey window. In both datasets, we assume Gaussian noise with a known power spectral density. In Figure 3, we plot these two simulated signals.

Figure 3. The correctly specified waveform (pink) and the misspecified waveform (grey) used in Section 3.2, plotted in the time domain.

It is worth pausing to distinguish this discussion of model misspecification from the similar-but-different topic of model selection. If we were discussing model selection, we would keep the data fixed and compare it two different signal models. However, misspecification occurs when the analyst has not conceived of the correct model to test. Therefore, since we are discussing misspecification, we keep the model fixed and consider two hypothetical datasets: one correctly specified and one misspecified.

In the right-hand panel of Figure 4a, we display the time-domain residuals obtained when we subtract the sine-Gaussian waveform model from the misspecified data. In the left-hand panel, we show the time-domain residuals obtained from subtracting the same waveform model from the correctly specified data. If the waveform is correctly specified, the residuals should be consistent with our Gaussian noise model as in the left panel. Although it is sometimes possible to see misspecification in the time-domain data (for example, when there is a short glitch), the misspecification may be more apparent in the frequency domain as is the case here.

Figure 4. Identifying a misspecified signal model. The left-hand column shows tests performed on data containing a signal consistent with the sine-Gaussian pulse model that we test against. The right-hand column shows the same tests performed on data containing a different signal.

In Figure 4b, we plot the amplitude spectral density of the residuals. Again, the left-hand panel shows the residuals for the correctly specified signal while the right-hand panel shows the residuals for the misspecified signal. While the correct waveform yields residuals consistent with the noise amplitude spectral density, the misspecified signal produces a peak in the amplitude spectral density that is conspicuously outside of the 90% range predicted by the model (shown in pink).

In order to assess the overall goodness of fit, it can be useful to plot the cumulative density function (CDF) of the residuals alongside the theoretical CDF predicted by the model. The CDF may be constructed from residuals in the time domain (to highlight transient phenomena), but more often, misspecification is most apparent using whitened frequency-domain residuals. An illustration of the CDF test is provided in Figure 4c, where we plot the CDF of the time-domain residuals. While there is only one realisation of the data (grey), we can generate arbitrarily many realisations of the theoretical CDF (pink). The thickness of the pink CDF shows the variability from 100 different realisations.(This is similar to, although not the same as, Gelman’s posterior predictive checks (e.g., Gelman et al. Reference Gelman, Carlin, Stern, Dunson, Vehtari and Rubin2013; Gelman & Rohilla Shalizi Reference Gelman and Rohilla Shalizi2010). Gelman’s checks involve drawing realisations from the histogram of the posterior probability distribution under the assumption of the model, and checking how probable it is for the model to produce realisations that are consistent with the observed data. This is something that we return to in Section 4. Here, we draw realisations from the noise model in order to establish the range over which it can credibly vary.)

To determine if the residuals agree with the model, one can employ any of the many established hypothesis-testing tools available to determine if measured samples are drawn from the theoretical distribution. For example, the Kolmogorov-Smirnov (KS) statistic,

(10) \begin{align} D_n = \max_f\Big| \text{CDF}(f|\text{measured}) - \text{CDF}(f|\text{predicted}) \Big| ,\end{align}

is the maximum difference between the measured CDF of the data and the predicted CDF given n frequency bins. (The Anderson-Darling test is also commonly used to determine if measured samples are drawn from a predicted distribution.) The $D_n$ statistic can be converted into a p-value with a look-up table that does not depend on the functional form of the CDF. The Kolmogorov-Smirnov p-values for the correctly specified data and misspecified dataset are provided in the legends of Figure 4c. The small p-value in the right panel suggests that our signal model is indeed misspecified.

Sometimes one may wish to inspect the data within a particular time or frequency interval. In this case, it can be enlightening to draw hundreds of distributions from the model, and count the number of times that the model CDF in this bin is above the CDF of the data. The fraction of draws above the CDF in this bin constitutes a p-value; if the data are correctly specified, the p-value will follow a uniform distribution. Thus, a p-value very close to 0 or 1 indicates that the data in this bin deviate significantly from the model. However, care must be taken if more than one p-value is calculated in the same set of data—while each p-value is uniformly distributed, the set of p-values from one dataset are correlated with each other.

Both of the tests described above assume that the predicted model is non-parametric. That is, the distribution of the residuals does not depend on any parameters. If the predicted distribution depends on one or more parameter $\theta$ , then we must rely on Monte Carlo methods to determine if the distributions agree. For example, we may calculate

(11) \begin{align} D'_n = \max_f\Big| \text{CDF}(f|\text{meas.}) - \text{CDF}(f|\widehat\theta,\text{pred.}) \Big| ,\end{align}

which is a KS-type statistic using the maximum-likelihood parameters $\widehat\theta$ . Since we use the data to estimate $\widehat\theta$ , $D'_n$ is not distributed according to the Kolmogorov distribution. (It is easier to get smaller differences between the measured and predicted CDFs when the theoretical CDF varies depending on the parameters.) However, we can still calculate a p-value by generating synthetic data to empirically estimate the distribution of $D'_n$ , which amounts to a generalisation of Lilliefors test (Lilliefors Reference Lilliefors1967). We demonstrate such a calculation below in Subsection 4.2 in the context of prior misspecification; see in particular Figure 7a.

3.3. Testing for a misspecified noise model

In order to observe noise misspecification and deviations from the Whittle likelihood, it is again useful to look at the distributions of residuals. By comparing the observed residuals to the theoretical likelihood distribution, it is possible to see, for example, if five-sigma deviations are more common than expected.

Recommended tests:

• • Create a histogram representing the probability density function of the data; or, alternatively, plot the CDF of $r,|\tilde{r}|$ , the whitened residuals. It is sometimes also useful to plot the distribution of the whitened residual power $|\tilde{r}|^2$ ; see Talbot & Thrane (Reference Talbot and Thrane2020), Chatziioannou et al. (Reference Chatziioannou, Haster, Littenberg, Farr, Ghonge, Millhouse, Clark and Cornish2019). Include residual curves for many posterior distribution draws of $\theta$ in order to show theoretical uncertainty.

• • For fixed models with no parameters, calculate a goodness-of-fit statistic using, for example, the KS test. When models have parameters, calculate the goodness of fit for the maximum-likelihood parameters.

• • Plot the difference in the cumulative density functions (empirical - expected) as a function of the expected CDF as in Talbot & Thrane (Reference Talbot and Thrane2020), Chatziioannou et al. (Reference Chatziioannou, Haster, Littenberg, Farr, Ghonge, Millhouse, Clark and Cornish2019). This is like a probability–probability (’PP’) plot, in which the fraction of repeated measurements is plotted against the confidence level at which the known truth value exists, for data model checking.

• • Bootstrapping methods—in which real data is used as a sampling distribution for synthetic data—can often be helpful for diagnosing noise misspecification. In gravitational-wave astronomy, new noise realisations can be created by bootstrapping data residuals (Cannon, Hanna, & Keppel Reference Cannon, Hanna and Keppel2013; Cannon, Hanna, & Peoples Reference Cannon, Hanna and Peoples2015; Ashton, Thrane, & Smith Reference Ashton, Thrane and Smith2019) or using time slides (Dal Canton et al. Reference Dal Canton2014; Usman et al. Reference Usman2016); both methods were integral to the first gravitational-wave detections (e.g., Abbott et al. Reference Abbott2016a). However, all bootstrap methods ultimately break down due to saturation effects (Wąs et al. 2010); it is impossible to simulate all possible noise realisations using a finite dataset. In astro-particle physics, using instruments such as Super-Kamiokande, sidereal time scrambling is used to estimate typical fluctuations in signal strength due to noise (Thrane et al. Reference Thrane2009).

Demonstration: Footnote g We demonstrate how to diagnose noise misspecification. The true noise is Gaussian with a mean $\mu=0$ and standard deviation $\sigma=1$ . The misspecified noise is distributed according to the Student’s t distribution with $\nu=5$ degrees of freedom. These parameters are chosen so that the noise profiles appear, at first glance, to be consistent with each other.

We display histograms of each noise dataset in Figure 5a. The region that the model predicts 90% of the data to lie within is also shown in these plots. The histograms both appear to largely lie within this 90% credible region, with the Student’s t-distributed data only visibly deviating in the tails.

Figure 5. Identifying a misspecified noise model. The left-hand column shows tests performed with a correctly specified Gaussian noise model while the right-hand column shows the same tests with the same Gaussian model, but performed against a misspecified Student’s-t distribution.

Next, we Fourier transform our datasets and compare the 90% range predicted by the noise model against histograms of the data in the frequency domain (Figure 5a). In the frequency domain, the misspecified data more clearly strays outside of the range predicted by the model. We then create a CDF of the frequency-domain data, comparing again to the 90% range predicted by the model (Figure 5b). We find a reasonable KS test p-value for the correctly specified data and an extreme p-value for the misspecified data, as expected. Finally, we compare the data to the model by plotting the data CDF against the difference between the data CDF and the model CDF (Figure 5c). In this final test, it is clear that the data is not well-represented by the model.

4.1. Priors with no hyper-parameters

In some cases, we can be very confident in our priors. For example, since there is no preferred direction in the Universe, the best prior for inclination angle $\iota$ (the angle between the orbital angular momentum and the line of sight) is uniform in $\cos\iota$ . In other situations, we are not confident in the form of the prior distribution. In these cases, it can be useful to formalise this theoretical uncertainty using a conditional prior

(12) \begin{align} \pi(\theta | \Lambda) .\end{align}

Here, $\Lambda$ is a ‘hyper-parameter’ we may vary to alter the shape of the prior for $\theta$ . In this Subsection, we focus on priors with no hyper-parameters while we cover parameterised priors in the next Subsection.

Recommended tests:

• • Make a CDF plot comparing the reconstructed distribution of $\theta$ to the expected distribution of $\theta$ . In order to obtain a reconstructed distribution, obtain posterior samples for N different events, each weighted with the population model to be tested. Draw one posterior sample from each of the N events to make a realisation of the reconstructed distribution. Do this many times to make many realisations. Plot CDFs of the many realisations alongside the population model being tested. If the data agrees with the model, the CDFs should overlap. If the two CDFs do not overlap, the model is a poor description of the data. It is worth noting that a model can pass this test while still being quite badly misspecified. Abbott et al. (Reference Abbott2021e) found evidence for negatively aligned spins in a population of compact-object binaries, but this was later shown to be a model-dependent feature (Roulet et al. Reference Roulet, Chia, Olsen, Dai, Venumadhav, Zackay and Zaldarriaga2021; Galaudage et al. Reference Galaudage, Talbot, Nagar, Jain, Thrane and Mandel2021). This test may be particularly unreliable if there is a sharp feature in the data that is not captured by the prior.

• • If a more quantitative test is desired, one may calculate the KS-statistic for each CDF. If the best-fitting reconstruction has a p-value below a certain threshold, for example p-value $\leq0.00005$ (Benjamin et al. Reference Benjamin2017), then the model is not a good fit to the data.

• • One may also identify problem areas of model under- or over-production by quantifying the discrepancies between the data draw CDFs and the model over the parameter space. This facilitates statements like ‘99% of the time, the model produces too many binary black hole events with $m_1>80\,{\rm M}_\odot$ ’.

Figure 6. Plots illustrating how model misspecification may manifest for an unparameterised population model, with a correctly specified dataset shown in the left-hand plots, and the misspecified case in the right-hand plots.

Demonstration. Footnote h We simulate data d consisting of some physical parameter x and noise n:

(13) \begin{align} d = x + n .\end{align}

The noise is Gaussian distributed with zero-mean and unit variance. The values of x are drawn from the true prior distribution: a Gaussian distribution with mean $\mu=10$ and width $\sigma=2$ . However, in order to demonstrate prior misspecification, we employ as our prior a Laplacian with the same mean and variance as the true distribution. The true prior distribution and the misspecified prior distribution are compared in Figure 6a.

We create a simulated dataset of $N=1000$ events. In each case, we construct Gaussian likelihood curves of mean $\mu_i = d_i$ . The posterior for each value of $x_i$ is proportional to the likelihood multiplied by the prior. We draw 50 posterior samples from each of these simulated posteriors in order to create 50 CDF curves; see Figure 6b. We make two versions of this plot: one with the correctly specified prior and one with the misspecified prior.

It is also sometimes useful to plot the data CDF minus the model CDF. We include an example of this plot in Figure 6c. In this case, we see that the misspecified prior yields disagreements in the CDF in the distribution tails; the difference between pink model and grey data is, in general, inconsistent with zero when the model CDF is ${\approx}0.1$ and ${\approx}0.9$ .

4.2. Hyper-parameterised priors

We are particularly concerned with population models (sometimes called ‘hierarchical models’), where the prior for the parameters $\theta$ is conditional on a set of hyper-parameters $\Lambda$ , describing the shape of the prior distribution:

(14) \begin{align} \pi(\theta | \Lambda) .\end{align}

While data models are built on assumptions about the nature of noise, population models are built on assumptions about astrophysics. For example, if we assume that the distribution of primary black hole mass $m_1$ follows a power-law distribution with spectral index $\alpha$ ,

(15) \begin{align} \pi(m_1 | \alpha) \propto m_1^\alpha ,\end{align}

then $\alpha\in\Lambda$ is a hyper-parameter for the distribution of $m_1 \in \theta$ .

Hyper-parameters are frequently a convenient way of describing systematic theoretical uncertainty. For example, consider a parameter x drawn from a Gaussian distribution ${\mathcal{N}}(x|\mu,\sigma)$ with mean $\mu$ and width $\sigma$ . If we do not know the precise value of $\mu$ or $\sigma$ , our misspecification tests should take this uncertainty into account. We repeat the tests from the previous Subsection for the case of a hyper-parameterised prior.

Here is a summary of the differences that arise from the addition of a hyper-parameter. One must generate random draws of the hyper-parameter $\Lambda$ . Then generate draws of the model CDF for each random draw of $\Lambda$ , and weight posterior samples using the population predictive distribution—the conditional prior, marginalised over uncertainty in the hyper-parameters:

(16) \begin{align} \text{ppd}(\theta) = \int d\Lambda\, \pi(\theta | \Lambda) p(\Lambda | d). \end{align}

Here $p(\Lambda | d)$ is the posterior for the hyper-parameters. Take care not to double-count; the posterior for $\Lambda$ used to reweight event k should not be informed by event k; see, for example, Essick et al. (Reference Essick, Farah, Galaudage, Talbot, Fishbach, Thrane and Holz2022).

Demonstration: Footnote i We assume a Gaussian prior for parameter x with uncertain $\mu$ and width $\sigma$ . We simulate $N=1000$ true events, $x_\mathrm{true}$ , from two populations (priors): one Gaussian-distributed and one Laplacian-distributed, both with the same means $\mu_\mathrm{G}$ , and with widths $\sigma_\mathrm{G}$ and $\sigma_\mathrm{L}$ , respectively. For each population, we calculate the maximum-likelihood detected value of each event, $x_\mathrm{meas}$ by offsetting it by a random number drawn from a Gaussian of the same width as the likelihood distribution, $\sigma_\mathrm{meas}$ . The mean average of these maximum-likelihood values is the maximum-likelihood estimate for the population prior mean, $\mu_\mathrm{E}$ , and the measured population prior width is $\sigma_\mathrm{E}=\sigma_\mathrm{meas} / \sqrt{N}$ . The estimated mean of the population posterior is described by

(17) \begin{equation} \mu_\mathrm{P} = \frac{\sigma_\mathrm{E}^{-2} \mu_\mathrm{E} + \sigma^{-2} \mu}{\left(\sigma_\mathrm{E}^{-2} + \sigma^{-2}\right)} = \frac{\mu_\mathrm{E}}{\sigma_\mathrm{E}^2 \left(\sigma_\mathrm{E}^{-2} + \sigma^{-2}\right)},\end{equation}

since $\mu = 0$ . The estimated width is

(18) \begin{equation} \sigma_\mathrm{P} = \left(\sigma_\mathrm{E}^{-2} + \sigma^{-2}\right)^{-\frac{1}{2}}.\end{equation}

The width of the posterior predictive distribution is

(19) \begin{equation} \sigma_\mathrm{pp} = \sqrt{\sigma_\mathrm{P}^{2} + \sigma^{2}},\end{equation}

the width of a cross product of the uncertainty distribution of width $\sigma$ and the posterior width $\sigma_\mathrm{P}$ . The population-weighted width for each individual event’s posterior is

(20) \begin{equation} \sigma_\mathrm{x} = \left(\sigma_\mathrm{E}^{-2} + \sigma_\mathrm{pp}^{-2}\right)^{-\frac{1}{2}}.\end{equation}

We draw $M=100$ samples from each of these posteriors of mean $x_\mathrm{meas}$ and width $\sigma_\mathrm{x}$ . In practice, we do this by drawing N samples M times from the population posterior. The CDFs of these draws can then be compared to the model CDF.

Since we have used the data to estimate the hyper-parameters $(\mu, \sigma)$ , we cannot use the standard KS test to ascertain the degree of misspecification. Instead, we must define our own KS-like test for the p-value using the distance between the model and data. We calculate 100 KS distances between the model and draws from the model. We can histogram this data to show the distribution of KS distances expected if the model is a good fit to the data; see the pink histogram in Figure 7a. We measure the KS distance between the average data realisation and the data-conditioned model, which is shown by a grey line in Figure 7a. The p-value equates to the fraction of the model KS distance distribution above the data KS distance. In Figure 7b, we plot the 90% credible bands of the parameterised model and compare against the 90% credible range of the data. Despite being conditioned on the same data, the misspecified model strays outside of the credible band of the data and achieves a p-value of only 0.01. Finally, in Figure, we show the same CDFs with the CDF of the model subtracted, and plot against the model CDF. This has the effect of emphasising the extent of the mismatch between the data and the misspecified model, while the well-specified model contains the median data draw over the entire range.

4.3. Models with more than one dimension

It is difficult to look for model misspecification in more than one or two dimensions. Tests like those described in this Article so far can be extended, but the tell-tale signs of misspecification (e.g., detectable structure in residuals) can be difficult to see in large-dimensional spaces. However, there are some tools available.

• • Make two-dimensional scatter plots. If the model does not have any parameters, it can be represented by contours while the data can be represented with two-dimensional error bars or credible intervals.

• • If the model is subject to theoretical uncertainty (i.e., it is described by hyper-parameters), it may be necessary to draw multiple contours. However, the plot can be difficult to read if there are more than two variables.

• • There is no standardised test for goodness of fit in two or more dimensions. Boutique tests designed to identify particular forms of misspecification can work well, but one must beware trial-factor penalties when designing the test after looking at the data.

Figure 7. Model misspecification for a parameterised population model, with the correctly specified case demonstrated in the left-hand side of the lower two plots, and the misspecified case demonstrated on the right.

5.1. Posterior stability

In the previous Sections we describe tests to find misspecification. However, misspecification is sometimes manifest from surprising results. One such example is a phenomenon we refer to as ‘posterior instability.’ Let us consider the posterior for some parameter (call it $\theta$ ) and imagine how this posterior changes as we accumulate data. On average, we expect the posterior to narrow as we include more information. It also tends to shift around, but only within the bounds of the previous credible regions. It would be surprising if the addition of data produced a posterior favouring $\theta=3$ when an earlier posterior (calculated with less data) disfavoured $\theta=3$ with high credibility.

In such cases we say that the posterior is not stable to the addition of data. This can be indicative of model misspecification. An example from gravitational-wave astronomy is the maximum black hole mass parameter, which proved to be unstable moving from GWTC-1 (Abbott et al. Reference Abbott2019b) to GWTC-2 (Abbott et al. Reference Abbott2020a), under the assumption that the primary black hole mass distribution is a power law with a sharp cut-off. When the mass model was improved to allow for additional features (deviations from a power law Talbot & Thrane Reference Talbot and Thrane2018), the maximum-mass parameter stabilised.Footnote j For another example of posterior instability from optical astronomy, see Liu et al. (Reference Liu, Gezari and Miller2018), Zhu & Thrane (Reference Zhu and Thrane2020).

5.2. Outliers

Outliers are symptom of misspecification closely related to posterior stability (Fishbach, Farr, & Holz Reference Fishbach, Farr and Holz2020; Essick et al. Reference Essick, Farah, Galaudage, Talbot, Fishbach, Thrane and Holz2022). An outlier is an event with a parameter value appearing inconsistent with the rest of the distribution—in the context of a particular population model. Let us imagine that the distribution of events in our dataset (characterised by parameter $\theta$ ) seem well-described by a normal prior with mean zero and unit variance: ${\mathcal{N}}(\theta|\mu=0, \sigma=1$ . If we subsequently observed an event with $\theta \gtrsim 10$ , this could indicate that our prior model (a normal distribution) is inadequate.

Commonly, population outliers are identified using a ‘leave-one-out’ analysis method that compares an inferred distribution with and without the potentially anomalous data point; see Fishbach et al. (Reference Fishbach, Farr and Holz2020). However, care must be taken to take into account trial factorsFootnote k since, in any catalogue, some event has to be the most extreme. The statistics required for a careful leave-one-out analysis are too complicated for us to summarise here. Instead we refer the reader to Essick et al. (Reference Essick, Farah, Galaudage, Talbot, Fishbach, Thrane and Holz2022), which describes a ‘coarse-graining’ method to identify outliers.