Estimating NPDE from Simulated Replicates

From Simulated Outcomes to NPDE

A predictive check can ask more than whether an observed value lies inside a simulated interval. It can ask where the observation falls in the model’s predictive distribution, then put that location on a common scale across observations. The normalized prediction distribution error (NPDE) does this using two operations: a rank-based normal transformation and a covariance adjustment.

The first operation expresses each observed outcome as a predictive rank. The second accounts for dependence among the normal scores of observations from the same subject. The calculations depend on having simulated joint replicate vectors under the fitted model: each replicate should preserve the model’s dependence structure across that subject’s observation times.

This construction follows the covariance transformation developed in the earlier material on deriving the NPDE covariance transformation. Here the emphasis is on estimating the ranks and covariance from a finite set of simulated replicates. For continuous outcomes, \(\Phi^{-1}\) is the standard normal quantile function. A rank close to 0 maps to a large negative score, and a rank close to 1 maps to a large positive score.

Estimate Ranks Without Infinite Scores

With finitely many simulations, the observed value can fall below every simulated value or above every simulated value. Using the raw empirical fractions \(0\) or \(1\) would produce \(\Phi^{-1}(0)=-\infty\) or \(\Phi^{-1}(1)=+\infty\). A finite-sample plotting position avoids these endpoints. One clear convention is to assign the observation the midrank it would have among the \(K\) simulated values and itself, then use the plotting position \((\text{midrank}-0.5)/(K+1)\).

Equivalently, let \(L_j\) be the number of simulated values strictly below \(y_j\), and let \(E_j\) be the number equal to \(y_j\). The midrank plotting position is

For a continuous outcome without ties, \(E_j=0\). If all \(K\) simulations are below the observation, then \(p_j=(K+0.5)/(K+1)\), which is less than 1; if they are all above it, \(p_j=0.5/(K+1)\), which is greater than 0. Thus every transformed score is finite. Ties can arise from rounded measurements, discrete outcomes, or repeated simulated values. Using the midrank gives equal outcomes equal ranks rather than imposing an arbitrary ordering among them.

The plotting-position choice should be stated and used consistently. Other rank corrections are possible, but different corrections can give different scores, especially when \(K\) is small or an observation is at an extreme. The choice does not change the fundamental sequence of operations: estimate a predictive rank, transform it to a normal score, and then account for covariance.

Estimate the Covariance from Simulated Vectors

The scores \(z_j\) for observed outcomes are marginal scores. They do not, by themselves, account for dependence between observation times. To estimate that dependence on the same rank-normal scale, transform each simulated replicate at each observation time. For component \(j\), rank \(Y_{1j},\ldots,Y_{Kj}\) among the \(K\) simulated values. If \(q_{kj}\) is the midrank of \(Y_{kj}\), use the finite positions

These positions lie strictly between 0 and 1. For each replicate \(k\), keep the scores across all \(m\) observation times together as one vector \(z_k^{\mathrm{sim}}\). Estimate \(V\) as their sample covariance matrix, using denominator \(K-1\). Its diagonal entries describe the estimated variances of the simulated rank-normal scores; its off-diagonal entries describe their estimated covariances. Because the original replicate vectors are joint simulations, this step retains the simulated dependence across times.

The covariance adjustment uses the symmetric inverse square root. If \(V=Q\Lambda Q^\mathsf{T}\), where \(Q\) contains orthonormal eigenvectors and \(\Lambda\) contains the positive eigenvalues, define

The matrix must be positive definite for this formula to be defined. In practice, a nearly zero eigenvalue is also a warning: the available simulations may not support a stable estimate of the covariance transformation.

Why the Transformation Removes Estimated Covariance

The key justification is a covariance identity. If \(z\) has covariance matrix \(V\), and \(W=V^{-1/2}\), then the covariance matrix after a linear transformation is \(WVW^\mathsf{T}\). Since the symmetric inverse square root is symmetric and \(V=Q\Lambda Q^\mathsf{T}\), direct substitution gives

Here \(I\) is the identity matrix. The calculation explains the purpose of whitening: under the covariance model used in the calculation, the transformed vector has identity covariance. With a finite number of replicates, \(V\) is estimated rather than known, so this result describes the transformation based on that estimate, not a guarantee that every finite-sample NPDE vector has exact identity covariance or a normal distribution.

Ties, Diagnostics, and Alternatives

A reliable NPDE calculation requires suitable simulations and a numerically usable covariance estimate. Before interpreting the result, check how the scores were constructed and whether the transformations behaved as intended.

• Check the simulated records. Each replicate must contain the corresponding simulated outcomes for the subject and observation times being analysed. Independently shuffling or separately simulating time points would remove the joint dependence that \(V\) is meant to estimate.
• Check rank extremes and ties. Record the rank convention, including its tie handling. Many observations at extreme ranks, or many ties caused by rounding, can make the result sensitive to the chosen plotting position.
• Check the covariance matrix. It should be symmetric and positive definite for the symmetric inverse square root. Very small eigenvalues can produce large transformed components; a non-positive eigenvalue prevents the ordinary inverse-square-root calculation.
• Check units and indexing. Rank calculations compare an observation with simulations for the same outcome and observation time. The rank-normal score and NPDE are dimensionless, even when the original outcome has units such as mg/L.

A common mistake is to stop after calculating \(z_j=\Phi^{-1}(p_j)\) and call that score an NPDE. That is only the marginal rank-normal score. When observations within a subject are correlated, the NPDE also requires the covariance adjustment \(V^{-1/2}z\). A different mistake is to rank each time point separately and discard which simulated values belonged to the same replicate; doing so loses the joint covariance information.

The symmetric inverse-square-root method is the standard choice here. A Cholesky factorization \(V=LL^\mathsf{T}\) offers another whitening transformation: \(L^{-1}z\) has identity covariance under the same covariance calculation. However, Cholesky whitening depends on the ordering of components and does not generally give the same transformed values as symmetric whitening. The two transformations serve the same covariance-removal purpose, but their component-wise results are not interchangeable.

NPDE also differ from visual predictive checks. A VPC summarizes observed and simulated quantiles across bins; NPDE produces transformed scores for individual observations and adjusts for within-subject covariance. Neither replaces the other. Rank and covariance estimates remain limited by the number and quality of simulations, by the chosen treatment of ties, and by whether the simulated records correctly represent the intended predictive distribution.