Identifiability, posterior predictive checks, and comparisons between ABC and SMM

Julia Fukuyama

Simulation-Based Inference · Lecture 15

Today: Wrapping up the basic ABC material, comparison to SMM

Just as with SMM:

  • We need informative summaries to use to match simulated and observed data, preferably sufficient statistics.
  • If you want to interpret estimated parameters, you need to check for identifiability
  • Once you have fit a model, you should check whether it reproduces the data. In the Bayesian context, this is a posterior predictive check.

In contrast to SMM:

  • ABC with a particular summary uses the whole distribution of the summary, not just the mean.

Recall: SMM and ABC

SMM: For each candidate parameter, average summaries across simulations. Find parameters whose averages match the observed summaries.

ABC: Draw parameters from a prior. Simulate a dataset at each draw. Retain parameters whose simulated summaries are close to the observed summaries.

Both compare summaries, but they make different comparisons.

1. Identifiability

A model is identifiable when different parameter values give different distributions of the observations.

A model is not identifiable if there are two parameter settings that give the same observed-data distribution.

Changing the fitting method cannot create information absent from the observations.

How to diagnose non-identifiability

  • In SMM, look for either a ridge of minimizing values, or separated minima with the same or similar objective values.
  • In ABC, look for a ridge of high posterior mass, or separated regions of substantial posterior probability.

These sorts of structures in the SMM objective or the ABC posterior give you candidates for non-identifiability, but they don’t prove that the model is not identifiable.

Investigating potential non-identifiability

Go back to the definition: the model is not identifiable if there are two parameter settings that give the same observed-data distribution.

Therefore:

  • Pick pairs of parameters along the ridge or pairs that correspond to different minima/posterior modes.
  • Simulate data from each.
  • If the distributions match, the model is not identifiable.
  • If the distributions don’t match, there is some true difference that your fitting procedure is not using.

Limitations of this method

  • A finite number of simulations cannot establish that two distributions are exactly the same.
  • The candidate parameter pairs are just guesses. Even if the model is exactly non-identifiable, picking pairs of parameter values from a ridge will almost always give slightly different distributions.

However, even showing that very different parameter values give very similar observed-data distributions tells us something valuable: If we really care about the parameters in the model, we need to either

  • Rely on very subtle differences in the data distributions, which means that estimates are sensitive to even small errors in modeling assumptions, or
  • Collect a different kind of data that distinguishes parameter values better.

Example: Insects caught in a trap

We count the insects caught each night.

We want to estimate the rate at which insects arrive in the area monitored by the trap, but some arrivals escape capture.

Let:

  • \(\lambda\) be the expected number of arrivals per night.
  • \(q\) be the probability of catching each arriving insect.

Could a small catch distinguish few arrivals from an ineffective trap?

The trap model

For each night \(i\):

\[ N_i\sim\operatorname{Poisson}(\lambda), \qquad Y_i\mid N_i\sim\operatorname{Binomial}(N_i,q). \]

and our prior is \[ \lambda &\sim \text{Unif}[0,40] \\ q &\sim \text{Unif}[0,1] \]

We observe the catches \(Y_1,\ldots,Y_n\). We do not observe the arrivals \(N_i\).

The retained parameter pairs

Draw \((\lambda,q)\) from the prior and simulate \(n\) nightly catches. Retain the parameter pair if the total simulated catch equals the observed total.

Left: 4,000 draws from the uniform prior. Right: the 3,702 of 2,000,000 proposals whose simulated 30-night total equals the observed total of 124. Both panels share axes.

This is a ridge with high posterior mass.

Step 1: Pick points along the ridge

The ridge suggests the model may not be identifiable. To check, pick parameter values far apart along it.

For \(\lambda\) near 5 and near 30, take the median \(q\) among retained pairs with \(\lambda\) within 10%.

The retained pairs from 30 nights, with the two chosen points marked.

Step 2: Simulate from each point and compare

Simulate many nights at each point and compare the two distributions of catches directly.

Catches on 100,000 simulated nights at each ridge point, with shared axes. Bins have width 1. Right: the two histograms overlaid.

The histograms match, up to simulation noise.

The same ambiguity in SMM

Using mean catch, the idealized SMM objective is

\[ Q(\lambda,q)=(\bar Y_{\mathrm{obs}}-\lambda q)^2. \]

Left: the exact SMM objective using the mean catch; the orange curve marks its minimizers, \lambda q = \bar y_{\mathrm{obs}} = 4.13. Right: the ABC posterior from 30 nights, with the same curve.

You could do the same thing with pairs of points along the minimizing ridge in SMM.

Would more data or more summaries help?

ABC posteriors using the first 30 nights and all 300 nights (exact match of the total catch, same prior). Right: the retained values of \lambda q.

More nights determine the product more precisely.

But the entire distribution of catches depends only on \(\lambda q\).

Even adding the variance, empty-night frequency, or other summaries cannot separate its factors.

Identifiability: What carries over?

As in SMM: We must ask whether the observations can distinguish the parameters we want to estimate.

With ABC: Inspect the joint posterior, not only separate parameter intervals. Different combinations may explain the same data.

More simulations improve computation. More observations sharpen what the data can identify. Additional information about capture or arrival is needed to separate \(\lambda\) and \(q\).

2. Model checking

Just as with SMM, a fitted model should reproduce relevant features of the observed data, including features we did not use for fitting.

Recall: Model checking with SMM

After fitting the model, we have a parameter estimate \(\hat\theta\).

To check the fit:

  1. Generate many new datasets from the model at \(\hat\theta\).
  2. Compare the simulated datasets with the observed dataset.
  3. Look for patterns in the observations that the simulations do not reproduce.

We can compare plots or choose a statistic \(T(X)\) and compare \(T(x_{\mathrm{obs}})\) with the simulated values.

Each simulated dataset should have the same sampling structure as the observed dataset.

Recall that you can have very precisely estimated parameters (a small bootstrap confidence interval), and still have a model that doesn’t fit well.

Bayesian version: Posterior predictive checks

With SMM, we simulated checking datasets at \(\hat\theta\).

With Bayesian inference, we have a posterior distribution \(p(\theta\mid x_{\mathrm{obs}})\).

To generate a new dataset, we first draw a parameter value from this posterior, then simulate data at that value:

\[ \theta^{(b)}\sim p(\theta\mid x_{\mathrm{obs}}), \qquad X_{\mathrm{rep}}^{(b)}\sim p(\,\cdot\mid\theta^{(b)}). \]

The distribution of these new datasets is the posterior predictive distribution.

The posterior predictive distribution

Let \(X_{\mathrm{rep}}\) be a new dataset generated independently of the observed dataset, conditional on \(\theta\).

Its distribution given the observations is

\[ p(x_{\mathrm{rep}}\mid x_{\mathrm{obs}}) = \int p(x_{\mathrm{rep}}\mid\theta)\, p(\theta\mid x_{\mathrm{obs}})\,d\theta. \]

This averages the model’s predictions over the posterior distribution of the parameters.

The simulation recipe on the previous slide performs this averaging without calculating the integral.

Posterior predictive checking with ABC

ABC gives us retained parameter draws approximating the posterior.

For each checking replicate:

  1. Choose one retained parameter draw.
  2. Simulate a fresh, complete dataset at that parameter value.
  3. Calculate the statistic \(T(X_{\mathrm{rep}})\).

Compare the observed statistic \(T(x_{\mathrm{obs}})\) with the distribution of replicated \(T(X_{\mathrm{rep}})\).

Does the fitted model commonly generate the feature we observed?

Example: Weather over several decades

We observe dry and rainy days during the same 30-day calendar period at one location, for 30 years.

We want to estimate:

  • \(p\): the overall probability of a dry day.
  • \(\rho\): how strongly weather persists from one day to the next.

Our initial model assumes both parameters are constant across years.

A model for the rain records

Let \(X_{ij}=1\) if day \(j\) of year \(i\) is dry and \(X_{ij}=0\) if it rains, for \(i,j=1,\ldots,30\).

Years are independent. Within year \(i\):

\[ X_{i1}\sim\operatorname{Bernoulli}(p), \qquad X_{ij}\mid X_{i,j-1}\sim\operatorname{Bernoulli}\big(\rho X_{i,j-1}+(1-\rho)p\big), \quad j=2,\ldots,30. \]

In words: each day copies yesterday’s weather with probability \(\rho\); otherwise it is a fresh draw, dry with probability \(p\).

What do the two parameters change?

\(p\) changes how often days are dry.

\(\rho\) changes how long wet and dry spells tend to last.

Even when we draw fresh weather, it may agree with yesterday.

Eight simulated years at each of four settings. Rows: p = 0.3 (top) and p = 0.7 (bottom). Columns: \rho = 0.1 (left) and \rho = 0.8 (right).

Each row is one year. Each column is one day.

Summaries for fitting

Use two summaries of \(X=(X_{ij})\):

  • The fraction of all days that are dry: \[ S_1(X)=\frac{1}{30\cdot 30}\sum_{i=1}^{30}\sum_{j=1}^{30}X_{ij}. \]
  • The fraction of within-year transitions that switch weather: \[ S_2(X)=\frac{1}{30\cdot 29}\sum_{i=1}^{30}\sum_{j=2}^{30}\mathbf 1\{X_{ij}\ne X_{i,j-1}\}. \]

The first informs dryness. Together, the two inform persistence.

ABC for the weather records

Priors: \(p\sim\operatorname{Uniform}(0,1)\) and \(\rho\sim\operatorname{Uniform}(0,1)\), independent.

For \(b=1,\ldots,B\):

  1. Draw \((p^{(b)},\rho^{(b)})\) from the prior.
  2. Simulate \(X^{(b)}\), 30 years of 30 days, from the rain model.
  3. Compute the standardized distance \[ d^{(b)}=\sqrt{\sum_{k=1}^{2}\left(\frac{S_k(X^{(b)})-S_k(x_{\mathrm{obs}})}{s_k}\right)^2}, \] where \(s_k\) is the standard deviation of \(S_k\) in a pilot run from the prior.

Retain the \((p^{(b)},\rho^{(b)})\) with the smallest 1% of distances.

The weather posterior

Left: joint ABC posterior keeping the closest 0.05% of 1,000,000 simulations (503 runs). Middle and right: marginal posteriors keeping the closest 1%, 0.2%, and 0.05%.

Observed: 52.1% of days dry, 28.1% of transitions switch. 90% intervals: \(p\) from 0.48 to 0.57, \(\rho\) from 0.39 to 0.49.

We obtain fairly narrow intervals. What assumptions do those intervals depend on?

Sampling from the posterior predictive distribution

With SMM, we simulated at the fitted parameter value.

With ABC, we have a collection of plausible parameter values.

For \(b=1,\ldots,R\):

  1. Choose a retained pair \((p^{(b)},\rho^{(b)})\).
  2. Simulate a fresh dataset \(X_{\mathrm{rep}}^{(b)}\) of 30 years from the rain model at that pair.

The resulting distribution of new datasets is an approximate posterior predictive distribution.

Observed weather and weather from the posterior predictive distribution

The observed records and the first three posterior predictive replicates, in generation order. Each row is a year, earliest at top; each column is a day.

Do the weather patterns drawn from the posterior predictive distribution look like the observed weather?

Turning a visual difference into a statistic

Let \(K_i\) be the number of weather switches in year \(i\):

\[ K_i=\sum_{j=2}^{30}\mathbf 1\{X_{ij}\ne X_{i,j-1}\}. \]

Our statistic is the least-squares slope of \(K_i\) against year \(i\):

\[ T(X)=\frac{\sum_{i=1}^{30}(i-\bar\imath)(K_i-\bar K)}{\sum_{i=1}^{30}(i-\bar\imath)^2}, \qquad \bar\imath=15.5,\quad \bar K=\frac{1}{30}\sum_{i=1}^{30}K_i. \]

A negative \(T\) means later years have fewer switches and longer spells.

Switches per year

Switches per year K_i for the observed records and the same three replicates, with least-squares lines; each panel title gives the slope T. All panels share axes.

Compare \(T(x_{\mathrm{obs}})\) with \(T(X_{\mathrm{rep}}^{(b)})\) for \(b=1,\ldots,R\).

Comparing the observed slope with posterior predictive slopes

T(X_{\mathrm{rep}}^{(b)}) for 4,000 posterior predictive datasets (histogram) and T(x_{\mathrm{obs}})=-0.36 (orange line).

Does the fitted model commonly generate a trend as strong as the observed one? None of the 4,000 replicates has \(T(X_{\mathrm{rep}}^{(b)})\le T(x_{\mathrm{obs}})\).

A mismatch here concerns the change across years, even if the pooled dry-day and switching fractions look reasonable.

Three things we want from our models that are not the same

  • Precise estimates of parameters (concentrated posterior or narrow bootstrap confidence interval)
  • Identifiable parameters
  • Model fits the data well (passes posterior predictive checks/replicates from the fitted look like the real data)

Non-identifiability will lead to at least some parameters not being estimated precisely, but precise parameter estimates and the model fitting the data well are completely different from each other and can vary independently.

Example: Passing posterior predictive checks does not establish identifiability

Return to the insect trap and its ABC posterior from 30 nights.

Check it as we checked the weather model: simulate fresh 30-night datasets, each at one retained \((\lambda,q)\) pair.

Compare features we did not use for fitting:

  • the variance of the nightly catches
  • the number of nights with no catch

Observed and replicated catches

Histograms of the observed catches (orange) and of the first 11 posterior predictive replicates, in generation order (blue). Each replicate is 30 nights simulated at one retained (\lambda,q) pair. All panels share axes.

Does the observed histogram stand out from the replicates?

The posterior predictive check

The variance of the catches and the number of empty nights in 8,000 posterior predictive datasets (histograms). Orange lines mark the observed values.

The observed variance (4.05) and number of empty nights (1) are typical of the replicated datasets.

The model passes this check.

3. SMM and ABC use summaries differently

We began with a procedural difference:

  • SMM matches expected summaries.
  • ABC retains individual simulations with matching summaries.

This means that ABC can use the entire distribution of a summary statistic rather than just the mean.

Silly example: A known mean and an unknown variance

Suppose

\[ X_1,\ldots,X_n\overset{\mathrm{iid}}{\sim}N(\mu_0,\sigma^2). \]

The population mean \(\mu_0\) is known. We want to estimate \(\sigma^2\).

For pedagogical reasons, assume we use only the sample mean \(S(X) = \frac{1}{n} \sum_{i=1}^n X_i\) as the summary.

Note: Yes, it is silly to try to use the sample mean to try to estimate the variance.

Does the distribution of the summary change with different parameter values?

Sample means of 500 simulated samples of size 25 at each standard deviation on a grid (points), with violins showing their distribution. Black bars: average of the simulated means at each \sigma. Dashed line: \mu_0 = 0.

The average of the simulated means is constant in \(\sigma\). Their variance grows with \(\sigma\).

Therefore: the mean of \(S(X)\) doesn’t contain any information about \(\theta\), but the distribution of the \(S(X)\) does.

SMM using the sample mean will not work

SMM compares the observed mean with the average of simulated sample means.

But

\[ E_{\sigma^2}[\bar X]=\mu_0 \qquad\text{for every }\sigma^2. \]

The idealized objective is therefore constant:

\[ Q(\sigma^2)=(\bar x_{\mathrm{obs}}-\mu_0)^2. \]

No variance fits better according to this objective.

Let’s try anyway

The SMM objective at 60 variances, estimated with 10 (gray) and 1,000 (blue) simulated samples per variance, and its exact value (orange).

With few simulations, the estimated objective has spurious minima. It flattens as we add simulations.

The noise also grows with the variance: those simulated means are more spread out.

ABC using the same sample mean

ABC retains a proposed variance if

\[ |\bar X_{\mathrm{sim}}-\bar x_{\mathrm{obs}}|<\epsilon. \]

The distributions of sample means from the previous figure. The orange area is the probability of landing within 0.05 of the observed mean, 0.4.

The probability of landing in that interval changes with the variance.

ABC therefore updates the prior even though the expected summary is constant.

The ABC posterior

ABC posteriors for \sigma^2 from 2,000,000 proposals at tolerances 0.05 and 0.02 (bars), the uniform prior (dashed), the exact posterior given the event that the sample mean lands within the tolerance (gray), and the variance used to simulate the observed sample (green).

The retained distribution reflects how often each variance generates a matching mean, as well as the prior.

The observed sample was simulated with \(\sigma^2 = 4\).

The full dataset would supply more: for this sample, \(\frac{1}{25}\sum_{i}(x_i-\mu_0)^2 = 3.64\).

Summing up

  • ABC and SMM both use summaries to estimate parameters, but they use the information in slightly different ways.
  • In each case, we have to check for non-identifiability and check for model fit.