Pooling Vector Estimates And Covariance Matrices

Why Pooling Vectors Requires More Than Pooling Standard Errors

Many analyses produce several estimates at once, such as a vector of regression coefficients. After multiple imputation, each completed data set gives one such vector and its estimated covariance matrix. Pooling each estimate separately can summarize its marginal uncertainty, but it does not retain how estimates vary together. That joint variation is recorded in a covariance matrix and matters whenever the estimates are considered jointly.

The matrix version of Rubin’s rules is a direct extension of scalar pooling. It averages the estimate vectors, averages their within-imputation covariance matrices, and adds a between-imputation covariance term with a finite-number correction. This tutorial focuses on those pooled quantities and what they represent; it does not develop a multivariate test.

Definition: Suppose there are \(m\) imputations, and each analysis estimates the same \(p\)-dimensional parameter vector. Let \(Q_i\) be the \(p\)-vector of estimates from imputation \(i\), and let \(U_i\) be its estimated \(p \times p\) within-imputation covariance matrix. Define the pooled vector estimate \(\bar{Q}\), average within-imputation covariance \(\bar{U}\), between-imputation covariance \(B\), and total covariance \(T\) by

\bar{Q}=\frac{1}{m}\sum_{i=1}^{m}Q_i,\qquad \bar{U}=\frac{1}{m}\sum_{i=1}^{m}U_i,

B=\frac{1}{m-1}\sum_{i=1}^{m}(Q_i-\bar{Q})(Q_i-\bar{Q})^{\mathsf T},\qquad T=\bar{U}+\left(1+\frac{1}{m}\right)B.

The transpose makes each outer product \((Q_i-\bar{Q})(Q_i-\bar{Q})^{\mathsf T}\) a \(p \times p\) matrix. Its diagonal entries are squared deviations for individual parameters; its off-diagonal entries are products of deviations and therefore describe how pairs of estimates vary across imputations. Both \(\bar{U}\) and \(B\) are covariance matrices, so their entries have units formed from the units of the relevant parameter estimates.

Assumptions and a Matrix Justification

Pooling is meaningful only if the \(Q_i\) vectors estimate the same \(p\) quantities in the same order and parameterization. Each \(U_i\) must be the estimated covariance matrix for that vector from the analysis of imputation \(i\), using the same model and covariance convention. For the usual inferential interpretation of Rubin’s rules, the imputations and analyses must also be appropriate to the missing-data problem: in particular, the imputation procedure should be proper, and the complete-data estimates should have a suitable approximately normal sampling distribution. The formula does not repair an incompatible or misspecified imputation or analysis model.

The terms in \(T\) have distinct roles. \(\bar{U}\) represents average uncertainty within completed-data analyses. \(B\) represents the additional variation in estimates across imputations, which reflects uncertainty about missing values. The factor \(1+1/m\) includes both that between-imputation variation and the extra uncertainty arising because only \(m\) imputations are used.

Justification: For any two coordinates \(j\) and \(k\), the \((j,k)\) entry of \(B\) is the sample covariance across imputations of the corresponding estimates. Thus, the same scalar pooling logic applies to every variance and covariance entry. The estimated total covariance adds the average within-imputation covariance to the between-imputation covariance, with the finite-\(m\) addition \(B/m\). Applying this entry by entry gives \(T=\bar{U}+B+B/m\), which is the stated matrix formula.

The matrices should be symmetric, up to small numerical rounding error. Since \(\bar{U}\) and \(B\) are covariance matrices, they are positive semidefinite in exact arithmetic; their weighted sum \(T\) is also positive semidefinite. If \(T\) is nonsingular, it can support calculations that require its inverse, but calculating or interpreting a particular joint test requires additional methods.

Worked Example: Pooling Two Estimates and Their Covariance

Suppose three hypothetical imputations produce two estimates: a treatment coefficient in log-odds units and an age coefficient in log-odds per year. The estimate vectors are \(Q_1=(1,2)^{\mathsf T}\), \(Q_2=(2,1)^{\mathsf T}\), and \(Q_3=(3,3)^{\mathsf T}\). Assume all analyses use this same order and parameterization, and their within-imputation covariance matrices are

U_1=\begin{pmatrix}0.16&0.04\\0.04&0.25\end{pmatrix},\quad U_2=\begin{pmatrix}0.25&0.05\\0.05&0.16\end{pmatrix},\quad U_3=\begin{pmatrix}0.09&0.03\\0.03&0.36\end{pmatrix}.

Pooled estimate. With \(m=3\), average each coordinate:

\bar{Q}= \begin{pmatrix}(1+2+3)/3\\(2+1+3)/3\end{pmatrix} = \begin{pmatrix}2\\2\end{pmatrix}.

Average within-imputation covariance. Average corresponding matrix entries:

\bar{U}= \begin{pmatrix} (0.16+0.25+0.09)/3&(0.04+0.05+0.03)/3\\ (0.04+0.05+0.03)/3&(0.25+0.16+0.36)/3 \end{pmatrix} = \begin{pmatrix}0.1667&0.0400\\0.0400&0.2567\end{pmatrix}.

Between-imputation covariance. The deviations from \(\bar{Q}\) are \((-1,0)^{\mathsf T}\), \((0,-1)^{\mathsf T}\), and \((1,1)^{\mathsf T}\). Their outer products sum to \(\begin{pmatrix}2&1\\1&2\end{pmatrix}\). Dividing by \(m-1=2\) gives

B=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}.

Total covariance. Here \(1+1/m=4/3\), so

T=\begin{pmatrix}0.1667&0.0400\\0.0400&0.2567\end{pmatrix} +\frac{4}{3}\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix} = \begin{pmatrix}1.500&0.7067\\0.7067&1.590\end{pmatrix}.

The diagonal entries are the total variances of the two pooled estimates. Their marginal standard errors are \(\sqrt{1.500}=1.225\) log-odds units and \(\sqrt{1.590}=1.261\) log-odds per year, respectively. The off-diagonal entry, \(0.7067\), is their total covariance, in the product of those units. It is not a variance or a standard error.

Reading and Checking the Covariance Matrix

A covariance matrix is not a table of independent standard errors. Its diagonal entries summarize uncertainty for individual coordinates; its off-diagonal entries summarize joint uncertainty. For coordinates \(j\) and \(k\), the total covariance can be inspected as

T_{jk}=\bar{U}_{jk}+\left(1+\frac{1}{m}\right)B_{jk}.

The within-imputation covariance \(\bar{U}_{jk}\) and between-imputation covariance \(B_{jk}\) can have different signs. A negative between-imputation covariance means that, across imputations, values of one estimate tend to be above their mean when values of the other tend to be below its mean. It does not by itself mean that the underlying parameters have a negative association in the population.

A useful descriptive summary of the pooled dependence is the correlation obtained from \(T\), when both diagonal entries are positive:

\operatorname{Corr}_{T}(j,k)=\frac{T_{jk}}{\sqrt{T_{jj}T_{kk}}}.

This correlation describes the dependence in the pooled estimation uncertainty, not necessarily the substantive relationship between the variables represented by the coefficients.

Worked Example: A Negative Between-Imputation Covariance

Four hypothetical imputations give estimates \((10,2)^{\mathsf T}\), \((12,1)^{\mathsf T}\), \((8,3)^{\mathsf T}\), and \((10,2)^{\mathsf T}\), in units appropriate to two distinct model coefficients. Suppose the average within-imputation covariance is

\bar{U}=\begin{pmatrix}4&1\\1&2\end{pmatrix}.

The pooled estimate is \(\bar{Q}=(10,2)^{\mathsf T}\). The deviations are \((0,0)^{\mathsf T}\), \((2,-1)^{\mathsf T}\), \((-2,1)^{\mathsf T}\), and \((0,0)^{\mathsf T}\). The sum of their outer products is \(\begin{pmatrix}8&-4\\-4&2\end{pmatrix}\); dividing by \(m-1=3\) gives

B=\begin{pmatrix}8/3&-4/3\\-4/3&2/3\end{pmatrix} \approx\begin{pmatrix}2.667&-1.333\\-1.333&0.6667\end{pmatrix}.

With \(m=4\), the correction factor is \(1.25\), and therefore

T=\begin{pmatrix}4&1\\1&2\end{pmatrix} +1.25\begin{pmatrix}2.667&-1.333\\-1.333&0.6667\end{pmatrix} \approx \begin{pmatrix}7.333&-0.6667\\-0.6667&2.833\end{pmatrix}.

The negative total covariance is valid: it records that uncertainty in the two estimates moves in opposite directions overall. The matrix is positive definite here, since its diagonal entries are positive and its determinant is approximately \(7.333(2.833)-(-0.6667)^2=20.33\), which is positive. A negative off-diagonal entry alone is not evidence of an invalid covariance matrix.

Changing Parameter Units

The pooling formulas behave coherently when the parameter vector is expressed in different linear units. If \(A\) is a fixed matrix and each estimate vector is re-expressed as \(Q_i^*=AQ_i\), the corresponding covariance matrices become \(U_i^*=AU_iA^{\mathsf T}\). Because averaging and matrix multiplication distribute over sums, the pooled quantities obey

\bar{Q}^*=A\bar{Q},\qquad \bar{U}^*=A\bar{U}A^{\mathsf T},\qquad B^*=ABA^{\mathsf T},\qquad T^*=ATA^{\mathsf T}.

Worked Example: Rescaling One Coordinate

Suppose a hypothetical analysis has \(\bar{Q}=(0.40,2.0)^{\mathsf T}\), \(m=4\), and

\bar{U}=\begin{pmatrix}0.04&0.01\\0.01&0.25\end{pmatrix},\qquad B=\begin{pmatrix}0.01&-0.005\\-0.005&0.04\end{pmatrix}.

First pool on the original scale:

T=\bar{U}+1.25B =\begin{pmatrix}0.0525&0.00375\\0.00375&0.3000\end{pmatrix}.

Now express the first coordinate in units ten times as large, leaving the second unchanged. This uses \(A=\begin{pmatrix}10&0\\0&1\end{pmatrix}\), so \(\bar{Q}^*=(4.0,2.0)^{\mathsf T}\). Transforming the total covariance gives

T^*=ATA^{\mathsf T} =\begin{pmatrix}5.250&0.03750\\0.03750&0.3000\end{pmatrix}.

The first variance has been multiplied by \(10^2\), and the covariance by \(10\), exactly as required by their units. The same result would be obtained by transforming \(\bar{U}\) and \(B\) separately and then applying the pooling formula. Mixing coefficient units or changing parameter order without applying the corresponding transformation would not be a valid comparison.

Diagnostics, Alternatives, and Limitations

Before using a pooled matrix, check that all estimates refer to the same parameters, in the same order, and on the same scale. Confirm that each \(U_i\) is a covariance matrix rather than a matrix of standard errors, and check that it is symmetric. Then verify that \(\bar{U}\), \(B\), and \(T\) are symmetric up to numerical rounding and have plausible diagonal entries. Materially negative variances or a substantially nonsymmetric result point to a calculation, coding, or model-output problem. A nearly singular matrix may also make later joint calculations unstable, even if its individual diagonal entries look reasonable.

The main alternatives answer different questions or discard information:

Pooling estimates but ignoring covariance. This can give the pooled point estimates, but it does not provide their joint uncertainty.
Pooling each coordinate separately. This can recover marginal variances, but it loses off-diagonal entries needed to describe dependence among estimates.
Using only \(\bar{U}\). This ignores uncertainty reflected by differences across imputations and generally understates total uncertainty when \(B\) is nonzero.
Analyzing one completed data set. This does not incorporate the between-imputation component and is not a substitute for pooling across imputations.

These matrix rules summarize estimates and covariance under Rubin’s pooling framework; they do not establish that the imputation model is suitable, guarantee reliable finite-sample approximations, or specify how to perform every joint inference procedure. When the number of imputations is small, the between-imputation covariance can itself be estimated imprecisely. A pooled covariance matrix should therefore be interpreted alongside the imputation design and diagnostics, not as a guarantee that all downstream inferences are reliable.

Common mistake: Treating the diagonal entries of \(T\) as standard errors, or setting its off-diagonal entries to zero because only marginal standard errors are reported. The diagonal entries are variances, so standard errors are their square roots. The off-diagonal entries are pooled covariances and should be retained when joint uncertainty matters.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

In one coordinate, estimates from three imputations are \(2, 4, 3\). Calculate the pooled estimate and that coordinate’s between-imputation variance.
Why is \(B\) divided by \(m-1\), rather than \(m\), in its definition?
If \(T_{11}=9\), what is the marginal standard error for the first coordinate, and what units does it have?
Can a negative off-diagonal entry of \(T\) be valid? Explain what it represents and name one check on the full matrix.
What must be true about the estimates before their vectors and covariance matrices can be pooled?

Solutions. (1) The mean is \(\bar{Q}=(2+4+3)/3=3\). The squared deviations are \(1,1,0\), so \(B=(1+1+0)/(3-1)=1\). (2) \(B\) is the sample covariance across the \(m\) imputation-specific estimates; dividing by \(m-1\) gives the usual sample covariance. (3) The standard error is \(\sqrt{9}=3\), in the units of the first parameter estimate. (4) Yes. It can indicate that the two estimates’ deviations tend to move in opposite directions. Check, among other things, that the whole matrix is symmetric and positive semidefinite, allowing for numerical rounding. (5) Every vector must estimate the same quantities, with the same ordering, parameterization, and units; each corresponding within-imputation covariance matrix must describe that vector’s uncertainty.