Constructing the Reference Predictive Distribution

What the Reference Predictive Distribution Represents

In reference-based imputation, a predictive distribution supplies plausible values for a participant’s outcome at scheduled visits. Its construction begins with a model for the reference arm: the fitted mean describes the expected outcome profile for a participant with specified covariates, and the fitted covariance describes the joint variation across visits around that profile. Combining the two gives a distribution for the participant’s outcome vector, not just a list of visit-specific averages.

This tutorial constructs that distribution as a plug-in multivariate normal model. The focus is on assembling and checking the reference-arm distribution. It does not condition on a participant’s observed post-baseline visits; that requires a different distributional calculation.

Definition: Let $Y_i=(Y_{i1},\ldots,Y_{iT})^\mathsf{T}$ be the $T$-visit outcome vector for participant $i$, and let $z_i$ be a $p$-dimensional vector of covariates, including an intercept if the model has one. A reference-arm mean model with common covariance is $$ Y_i=B^\mathsf{T}z_i+\varepsilon_i,\qquad \varepsilon_i\sim N_T(0,\Sigma_R). $$ Here $B$ is a $p\times T$ coefficient matrix and $\Sigma_R$ is a $T\times T$ covariance matrix for outcome variation within a participant in the reference arm. For a participant with covariates $z$, the plug-in reference predictive distribution is $$ Y_{\mathrm{new}}\mid z,\widehat B,\widehat\Sigma_R \ \sim\ N_T\!\left(\widehat B^\mathsf{T}z,\widehat\Sigma_R\right). $$

The first argument of this multivariate normal distribution is a vector of visit-specific fitted means; the second is a covariance matrix. In particular, the diagonal entries of $\widehat\Sigma_R$ are visit-specific variances, and its off-diagonal entries are covariances between visits. The distribution is participant-specific through $z$, even though this model uses a covariance matrix shared by reference-arm participants.

Assumptions and Construction

The model makes substantive and statistical commitments. Its mean specification must adequately describe how the expected reference-arm outcome profile varies with the included covariates. Conditional on those covariates, participants’ outcome vectors are assumed independent, with a common covariance matrix $\Sigma_R$. The multivariate normal assumption describes the joint residual distribution; it supports the stated normal predictive distribution and can be assessed, though not proved, from the data. The covariance matrix must be symmetric and positive definite for a nondegenerate multivariate normal distribution.

In practice, $\widehat B$ and $\widehat\Sigma_R$ are estimated using reference-arm data. One simple estimation setting has $n_R$ participants with complete outcome vectors and a common design matrix $Z$, whose rows are $z_i^\mathsf{T}$. With $Y$ the $n_R\times T$ outcome matrix, ordinary least squares gives

\widehat B^\mathsf{T}=(Z^\mathsf{T}Z)^{-1}Z^\mathsf{T}Y,\qquad E=Y-Z\widehat B^\mathsf{T},\qquad \widehat\Sigma_R=\frac{E^\mathsf{T}E}{n_R-p}.

This covariance estimator is unbiased under the fixed-design multivariate linear model with independent, mean-zero errors and common covariance $\Sigma_R$. The denominator is the residual degrees of freedom after fitting $p$ regression coefficients for each visit. If the reference-arm data are incomplete, the simple residual formula need not apply directly; the parameter estimates must instead come from an estimation procedure appropriate to the observed data and its assumptions.

Once the estimates are available, construction is direct: evaluate the fitted mean profile at the participant’s covariates, then pair it with the estimated reference-arm covariance matrix. Keeping the two components separate helps avoid a common error: covariates determine the mean in this model, but do not automatically change the covariance.

Fit the reference mean model.
Estimate the visit-specific regression coefficients using reference-arm data and the chosen covariates.

Estimate the residual covariance.
Use residual outcome vectors to estimate the variances and cross-visit covariances, checking that the result is a valid covariance matrix.

Evaluate at the participant’s covariates.
Compute $\widehat B^\mathsf{T}z$ and combine it with $\widehat\Sigma_R$ to specify the full plug-in predictive distribution.

Worked Calculations

Worked Example: Predicting Two Reference-Arm Visits

Suppose a hypothetical reference-arm analysis has estimated a two-visit mean of 10 and 8 units for a participant with specified covariates, and estimated covariance matrix $$ \widehat\Sigma_R= \begin{pmatrix} 4&2\\ 2&9 \end{pmatrix}. $$ Assume the fitted model’s mean and covariance assumptions are being used for this prediction. Construct the distribution and interpret its components.

The predictive mean vector is $(10,8)^\mathsf{T}$. The two marginal variances are 4 and 9 units squared, so the marginal standard deviations are $\sqrt{4}=2$ units and $\sqrt{9}=3$ units. The covariance between the visits is 2 units squared. The correlation is $$ \widehat\rho_{12}=\frac{2}{\sqrt{4}\sqrt{9}}=\frac{2}{6}=0.3333. $$ Thus the plug-in distribution is $$ Y_{\mathrm{new}}\mid z,\widehat B,\widehat\Sigma_R \sim N_2\!\left( \begin{pmatrix}10\\8\end{pmatrix}, \begin{pmatrix}4&2\\2&9\end{pmatrix} \right). $$

This specifies a joint distribution: the visit outcomes vary together, rather than being generated as independent values. The positive covariance means that, under this model, deviations above the fitted mean at one visit tend to accompany deviations above the fitted mean at the other. It does not mean that every participant follows that pattern.

Worked Example: Estimating Covariance From Residual Vectors

Consider hypothetical complete two-visit reference-arm data from three participants after fitting an intercept-only mean at each visit. Suppose the residual vectors are $(1,2)^\mathsf{T}$, $(-1,0)^\mathsf{T}$, and $(0,-2)^\mathsf{T}$, in outcome units. The residuals sum to zero at each visit, as expected for this intercept-only fit. Estimate the common residual covariance.

There are $n_R=3$ participants and $p=1$ fitted mean coefficient per visit, so the denominator is $n_R-p=2$. The residual cross-product matrix is $$ E^\mathsf{T}E= \begin{pmatrix} 1^2+(-1)^2+0^2 & 1(2)+(-1)(0)+0(-2)\\ 1(2)+(-1)(0)+0(-2) & 2^2+0^2+(-2)^2 \end{pmatrix} = \begin{pmatrix} 2&2\\ 2&8 \end{pmatrix}. $$ Therefore, $$ \widehat\Sigma_R=\frac{1}{2} \begin{pmatrix} 2&2\\ 2&8 \end{pmatrix} = \begin{pmatrix} 1&1\\ 1&4 \end{pmatrix}. $$

The estimated standard deviations are 1 and 2 outcome units, and the estimated correlation is $1/(1\cdot2)=0.5$. The matrix is positive definite: its leading diagonal entry is 1, and its determinant is $1(4)-1^2=3>0$. Thus it is a mathematically valid covariance matrix for a nondegenerate bivariate normal model. With a fitted mean vector of $(10,8)^\mathsf{T}$ for a specified participant, this covariance matrix completes the plug-in predictive distribution.

Worked Example: Incorporating a Participant Covariate

Suppose a hypothetical two-visit reference-arm mean model uses an intercept and age centered at 50 years. The estimated coefficient matrix is $$ \widehat B= \begin{pmatrix} 12&15\\ 0.5&0.4 \end{pmatrix}, $$ where the first row gives the intercepts and the second row gives the age slopes in outcome units per year. For a participant aged 54 years, $z=(1,4)^\mathsf{T}$. Suppose the estimated covariance is $$ \widehat\Sigma_R= \begin{pmatrix} 9&6\\ 6&16 \end{pmatrix}. $$ Find the predictive mean and the distribution of the difference between visit 2 and visit 1.

The participant’s fitted mean vector is $$ \widehat B^\mathsf{T}z= \begin{pmatrix} 12&0.5\\ 15&0.4 \end{pmatrix} \begin{pmatrix}1\\4\end{pmatrix} = \begin{pmatrix} 12+0.5(4)\\ 15+0.4(4) \end{pmatrix} = \begin{pmatrix}14\\16.6\end{pmatrix}. $$ The predicted mean is 14 outcome units at visit 1 and 16.6 outcome units at visit 2. The marginal standard deviations are $\sqrt{9}=3$ and $\sqrt{16}=4$ outcome units.

Let $D=Y_{\mathrm{new},2}-Y_{\mathrm{new},1}$. Its fitted mean is $16.6-14=2.6$ outcome units. For a contrast $a^\mathsf{T}Y_{\mathrm{new}}$, its variance is $a^\mathsf{T}\widehat\Sigma_Ra$; here $a=(-1,1)^\mathsf{T}$. Hence $$ \operatorname{Var}(D)=9+16-2(6)=13,\qquad SD(D)=\sqrt{13}=3.606\text{ outcome units}. $$ Since a linear combination of a multivariate normal vector is normal, the plug-in distribution of $D$ is $N(2.6,13)$, with the second argument denoting variance. This calculation uses the cross-visit covariance; adding the two variances without subtracting twice the covariance would give an incorrect variance for the difference.

Why the Joint Distribution and Its Checks Matter

The covariance structure affects more than the marginal uncertainty at each visit. It determines the uncertainty of changes, averages, and other combinations of outcomes. For any fixed vector $a$, the model implies $$ a^\mathsf{T}Y_{\mathrm{new}}\mid z,\widehat B,\widehat\Sigma_R \sim N\!\left(a^\mathsf{T}\widehat B^\mathsf{T}z,\ a^\mathsf{T}\widehat\Sigma_Ra\right). $$ For example, a mean across visits uses a vector $a$ with entries $1/T$; its variance depends on every entry of the covariance matrix, not only the diagonal. Ignoring cross-visit covariances can therefore distort uncertainty for summaries of a participant’s trajectory.

Useful diagnostics include checking that the estimated covariance is symmetric, has nonnegative variances, and is positive definite (or positive semidefinite in a model permitting degeneracy). Inspect residual means and residual patterns by visit and by important covariates. Systematic residual trends suggest that the mean model may be inadequate; markedly different residual spreads across visits or participant groups may challenge a common covariance assumption. Strongly non-normal residual patterns or influential participants can also make the multivariate normal approximation questionable, especially with limited reference-arm data.

A covariance estimate that is nearly singular deserves attention: it can make some simulated outcome combinations unstable and may signal too many covariance parameters relative to the available data. Check the number of reference participants, the rank of the design matrix, and whether the fitted model is trying to estimate more structure than the data can support. A numerical algorithm returning a matrix is not, by itself, evidence that the underlying assumptions are reasonable.

Important distinction: The plug-in distribution uses estimated parameters as though they were fixed. Its covariance $\widehat\Sigma_R$ represents participant-level outcome variation conditional on the fitted parameters; it does not, by itself, include uncertainty in $\widehat B$ or $\widehat\Sigma_R$. A predictive procedure that draws parameters or otherwise propagates their estimation uncertainty is a different construction and may produce wider uncertainty.

The marginal distribution constructed here is also distinct from a distribution conditional on a participant’s observed visits. The marginal model describes the full outcome vector for a participant with given covariates. Conditioning on some components changes the distribution for the others and requires additional calculations; substituting the marginal mean and covariance without conditioning would generally be wrong.

Compared with fitting a separate mean and variance at every visit, a multivariate model preserves cross-visit dependence and can represent covariate effects consistently across the outcome profile. Compared with a model that pools reference and non-reference outcomes, a reference-arm model uses the reference group to define its mean and covariance; pooling can change both. The choice of model is an assumption about the distribution to be used, not a fact established merely because the fitted matrix is numerically well behaved.

Common Mistakes and Limitations

Treating visit outcomes as independent. Using only the diagonal of $\widehat\Sigma_R$ discards the modeled cross-visit association. Use the full covariance matrix whenever a joint predictive distribution is required.
Confusing covariance with correlation. A covariance is measured in squared outcome units, while a correlation is unitless. Divide covariance by the product of standard deviations to obtain correlation.
Using a population-average mean for every participant. If the model includes covariates, evaluate its fitted mean at the participant’s covariate vector. Do not omit covariate contributions.
Calling a plug-in distribution fully uncertainty-adjusted. The plug-in distribution includes the estimated residual outcome variation but treats fitted parameters as fixed. Parameter uncertainty needs explicit propagation.
Assuming a valid-looking covariance proves a good model. Positive definiteness is a mathematical check, not evidence that the mean model, normality, common covariance, or reference-arm choice is appropriate.

The method is limited by the reference-arm sample size, covariate overlap, the chosen mean and covariance structures, and the adequacy of the distributional assumptions. Predictions for covariate profiles poorly represented in the reference arm rely heavily on model extrapolation. These limitations should be reported alongside any use of the resulting predictive distribution.

Check Your Understanding

Answer each question, then compare with the solutions below.

A two-visit model has fitted mean $(5,7)^\mathsf{T}$ and covariance matrix $\begin{pmatrix}4&-1\\-1&9\end{pmatrix}$. What are the marginal standard deviations and the correlation?
For a two-visit outcome, why is the variance of $Y_2-Y_1$ not generally the sum of the two marginal variances?
What assumption does a common $\Sigma_R$ make about participants with different covariates?
Does the plug-in predictive distribution automatically include uncertainty in the estimated regression coefficients? Explain.

Solutions. (1) The standard deviations are $\sqrt{4}=2$ and $\sqrt{9}=3$ outcome units. The correlation is $-1/(2\cdot3)=-0.1667$. (2) The variance is $\operatorname{Var}(Y_2-Y_1)=\operatorname{Var}(Y_2)+\operatorname{Var}(Y_1)-2\operatorname{Cov}(Y_1,Y_2)$; the covariance term is zero only if the visits are uncorrelated. (3) Conditional on the modeled covariates, participants share the same within-participant residual covariance, even though their fitted means may differ. (4) No. The plug-in calculation treats estimated coefficients as fixed; parameter uncertainty must be propagated separately.