Checking Candidate Families With Predictive Replication

Why Replicate Responses?

A fitted model summarizes how the response should behave at the observed predictor values. A useful check is to generate new responses from that fitted model, keeping those predictor values fixed, and ask whether the new data resemble what was observed. This is predictive replication: it assesses whether a candidate distribution can reproduce features such as exact zeros, variability, asymmetry, and unusually large outcomes.

The point is not to ask whether a simulated data set matches the observed data value by value. Replications vary by design. Instead, compare summaries calculated in the same way for the observed responses and for many replicated data sets. A single replication can be revealing, but a collection of replications shows how much variation the fitted model itself predicts.

For a GLM with fitted means \(\hat\mu_i\) and fitted dispersion \(\hat\phi\), the responses are generated from the specified candidate family at those values. The models considered here use \(\operatorname{Var}(Y_i\mid x_i)=\phi V(\mu_i)\): \(V(\mu)=1\) for Gaussian, \(V(\mu)=\mu^2\) for gamma, \(V(\mu)=\mu^3\) for inverse Gaussian, and \(V(\mu)=\mu^p\) for Tweedie. For a semicontinuous Tweedie model, \(1<p<2\). These variance functions describe conditional variability; replication checks ask whether the full generated distributions also reproduce important features of the data.

What Gets Simulated?

Each candidate model uses the same observed predictor values and its own fitted parameters. The response-generation step differs by family. In the following parameterizations, the fitted means must be positive for gamma, inverse-Gaussian, and semicontinuous Tweedie responses.

For the inverse-Gaussian parameterization, the density at \(y>0\) is

For the Tweedie model with \(1<p<2\), one way to generate a response is to draw \(N_i\sim\operatorname{Poisson}(\lambda_i)\), where \(\lambda_i=\mu_i^{2-p}/\{\phi(2-p)\}\), and sum \(N_i\) independent gamma jumps. Each jump has shape \((2-p)/(p-1)\) and scale \(\phi(p-1)\mu_i^{p-1}\). If \(N_i=0\), the sum is exactly zero. In particular,

Summaries to Compare

Choose summaries before inspecting the replications, and calculate each summary identically for the observed responses and every replicated data set. The summaries should reflect the feature being checked. For example, use the observed zero count when zeros matter, a variance or standard deviation for spread, a third-moment skewness for asymmetry, and a maximum or high quantile for the upper tail. When fitted means differ substantially across observations, also compare summaries within meaningful fitted-mean groups: a single overall variance can hide differences in expected variability across the design.

For a data vector \(y_1,\ldots,y_n\), common descriptive quantities include the exact-zero fraction, empirical variance, standardized third central moment, and maximum:

The displayed variance uses divisor \(n\), as do the central moments in this descriptive skewness. Other conventions are possible, but observed and replicated summaries must use the same convention. These raw summaries are easiest to interpret when observations have similar fitted means. With varying means, a discrepancy might reflect where observations sit in the design as well as the family’s shape; comparisons by fitted-mean group help make that distinction visible.

Proof. Define \(A_i=I(Y_i^{\mathrm{rep}}=0)\). Each \(A_i\) is Bernoulli with probability \(q_i\), so \(E[A_i]=q_i\) and \(\operatorname{Var}(A_i)=q_i(1-q_i)\). Linearity of expectation gives \(E[Z^{\mathrm{rep}}]=\sum_iq_i\). Conditional independence of the replicated responses implies independence of the indicators, so all their pairwise covariances are zero. Therefore \(\operatorname{Var}(Z^{\mathrm{rep}})=\sum_i\operatorname{Var}(A_i)=\sum_iq_i(1-q_i)\). When all probabilities equal \(q\), the sum of the independent, identically distributed Bernoulli indicators is binomial. \(\square\)

Using a Replicated Distribution, Not Just One Sample

Generate \(R\) independent replications, calculate a chosen summary \(T\) for each, and compare the observed value \(T_{\mathrm{obs}}\) with the resulting distribution of \(T^{\mathrm{rep}}\). A histogram or an interval covering the central bulk of the replicated summaries gives a practical reference. For example, if the observed zero count is far above nearly all replicated zero counts, the fitted model rarely produces as many zeros as were observed. If the observed maximum is larger than almost all replicated maxima, the model rarely produces a tail as extreme as the observed one.

The reference is conditional on the fitted parameters: it answers what the fitted model predicts when those parameter estimates are held fixed. If \(T\) has finite variance under that replication distribution, the simulation average \(\bar T_R\) estimates its conditional expectation. Independent replications give

These statements follow from linearity of expectation and the zero covariances between independent replicated summaries. The second equation also shows why more replications reduce Monte Carlo error, though they do not remove uncertainty arising because the model parameters were estimated from data.

Limits and a Common Mistake

A replication check asks whether a fitted model can plausibly reproduce selected features of the observed responses under the chosen predictors and fitted parameters. It does not identify the cause of a mismatch. A discrepancy might involve the response family, the mean specification, or another aspect of the data-generating process. Nor does a good match for one summary establish that all aspects of the distribution are adequate. Check several features that matter for the application, and interpret them together.

• Mistake: treating a plug-in replication as a new observed sample with known parameters. The replication holds fitted values fixed, so it does not include uncertainty from estimating those values.
• Mistake: declaring a model adequate because one simulated sample looks similar. Replications naturally vary; compare the observed summaries with their distributions across many generated samples.
• Mistake: expecting a continuous positive family to generate exact zeros. Gamma and inverse-Gaussian responses are strictly positive. Their zero probability is zero, regardless of how small a positive response may be.
• Mistake: comparing only a pooled spread when fitted means vary. Differences across the predictor design can affect the pooled summary; compare within fitted-mean groups as well.

The practical value of predictive replication is that it turns an abstract family choice into observable consequences. Candidate families make different predictions about zeros, asymmetry, variability, and extremes. Generating responses from each fitted candidate makes those predictions directly comparable with the data at hand.