Unbiasedness and Sampling Variation of Multivariate OLS

Why Coefficients Vary Across Samples

A multivariate regression records several responses on each case and uses the same design matrix to model them. Even when the regression mean is correctly specified, a new sample of cases will generally produce different estimated coefficients. Two features determine that sampling variation: the arrangement of the predictors and the covariance among responses within a case.

The response covariance can make estimates for different outcomes move together, while the design can make estimates for different predictors move together or in opposite directions. Neither feature by itself creates bias. The calculations below separate the expected coefficient matrix from its sampling covariance and make clear which assumptions support each conclusion.

Model and Notation

Definition: Let $Y$ be an $n\times q$ response matrix, $X$ a fixed $n\times p$ design matrix of full column rank, and $B$ a $p\times q$ coefficient matrix. The multivariate linear model is $Y=XB+E$, where the $i$th row $E_i$ of $E$ has mean zero and covariance matrix $\Sigma$, and the error rows are independent. Assume all error components have finite second moments. The ordinary least-squares estimator is $\hat B=(X^\top X)^{-1}X^\top Y$.

The covariance matrix $\Sigma$ is $q\times q$, with entry $\Sigma_{jk}=\operatorname{Cov}(E_{ij},E_{ik})$. It describes dependence among the $q$ responses within a case; it is the same for every row. Independence of error rows means errors from distinct cases have covariance zero. Normality is not required for the expectation or covariance results below.

Write $C=(X^\top X)^{-1}$, and let $C_{rs}$ denote its $(r,s)$ entry. For clarity, $\hat B_{rj}$ is the estimate for predictor $r$ and outcome $j$. The vectorization convention used here stacks the columns of a matrix: $\operatorname{vec}(\hat B)$ lists the coefficients for outcome 1, then outcome 2, and so on.

Mean and Covariance of the OLS Coefficients

Theorem: Under the model conditions above, $E[\hat B]=B$. For any predictors $r,s$ and outcomes $j,k$, $$ \operatorname{Cov}(\hat B_{rj},\hat B_{sk})=C_{rs}\Sigma_{jk}. $$ Equivalently, $$ \operatorname{Cov}\!\left(\operatorname{vec}(\hat B)\right)=\Sigma\otimes C. $$ These statements require zero-mean errors and the stated second-moment conditions; they do not require normally distributed errors.

Proof. Substitute $Y=XB+E$ into the least-squares formula:

\hat B=(X^\top X)^{-1}X^\top(XB+E) =B+(X^\top X)^{-1}X^\top E.

Since $X$ is fixed, taking expectations and using $E[E]=0$ gives $E[\hat B]=B$. For the covariance, set $A=(X^\top X)^{-1}X^\top$. The $(r,i)$ entry of $A$ gives the weight on case $i$ in coefficient row $r$, so

\hat B_{rj}-B_{rj}=\sum_{i=1}^n A_{ri}E_{ij}.

Expanding the covariance of two such sums gives

\operatorname{Cov}(\hat B_{rj},\hat B_{sk}) =\sum_{i=1}^n\sum_{l=1}^n A_{ri}A_{sl} \operatorname{Cov}(E_{ij},E_{lk}).

When $i\ne l$, independence of the error rows and their zero means imply $\operatorname{Cov}(E_{ij},E_{lk})=0$. When $i=l$, that covariance is $\Sigma_{jk}$. Thus the double sum reduces to

\Sigma_{jk}\sum_{i=1}^n A_{ri}A_{si} =\Sigma_{jk}(AA^\top)_{rs} =\Sigma_{jk}C_{rs},

because $AA^\top=(X^\top X)^{-1}X^\top X(X^\top X)^{-1}=C$. Under column-stacking vectorization, this entry-by-entry formula is exactly $\operatorname{Cov}(\operatorname{vec}(\hat B))=\Sigma\otimes C$. $\square$

In particular, the variance of one estimated coefficient and the covariance of two estimates are

\operatorname{Var}(\hat B_{rj})=C_{rr}\Sigma_{jj}, \qquad \operatorname{Cov}(\hat B_{rj},\hat B_{sk})=C_{rs}\Sigma_{jk}.

The first equation factors uncertainty into a predictor-design contribution and an outcome-variation contribution. The second shows that cross-outcome dependence and dependence induced by the design both matter when comparing two coefficients.

Worked Calculations

Worked Example: A Two-Predictor, Two-Outcome Model

Suppose the fixed design has $X^\top X=\begin{pmatrix}4&6\\6&14\end{pmatrix}$, and the independent error rows have a bivariate normal distribution with covariance $\Sigma=\begin{pmatrix}4&1.2\\1.2&9\end{pmatrix}$. Calculate the covariance matrix of the four OLS coefficients, in column-stacking order.

First invert the design cross-product matrix. Its determinant is $4(14)-6(6)=20$, so

C=(X^\top X)^{-1} =\frac{1}{20}\begin{pmatrix}14&-6\\-6&4\end{pmatrix} =\begin{pmatrix}0.7&-0.3\\-0.3&0.2\end{pmatrix}.

The coefficient vector is $(\hat B_{11},\hat B_{21},\hat B_{12},\hat B_{22})^\top$. Its covariance matrix is $\Sigma\otimes C$:

\begin{pmatrix} 2.8&-1.2&0.84&-0.36\\ -1.2&0.8&-0.36&0.24\\ 0.84&-0.36&6.3&-2.7\\ -0.36&0.24&-2.7&1.8 \end{pmatrix}.

For example, $\operatorname{Var}(\hat B_{11})=0.7(4)=2.8$, and $\operatorname{Var}(\hat B_{22})=0.2(9)=1.8$. The covariance between the two outcome-specific estimates for predictor 1 is $0.7(1.2)=0.84$. In contrast, the covariance between the two predictor estimates for outcome 1 is $(-0.3)(4)=-1.2$. The negative sign comes from the off-diagonal entry of $C$, not from a negative outcome covariance.

Worked Example: Intercept-Only Estimates for Correlated Outcomes

For $n=4$ cases, fit an intercept-only model to two outcomes. Suppose the error rows are independent bivariate normal random vectors with covariance $\Sigma=\begin{pmatrix}1&0.6\\0.6&4\end{pmatrix}$. Find the covariance matrix and standard errors of the two estimated means.

Here $X$ is a column of four ones, so $X^\top X=4$ and $C=1/4$. The two entries of $\hat B$ are the sample means for outcomes 1 and 2. Applying the coefficient covariance formula gives

\operatorname{Cov}(\operatorname{vec}(\hat B)) =\frac14\Sigma =\begin{pmatrix}0.25&0.15\\0.15&1\end{pmatrix}.

The standard errors are the square roots of the diagonal entries: $\sqrt{0.25}=0.5$ and $\sqrt{1}=1$. The covariance of the two estimated means is $0.6/4=0.15$. Their correlation is

\frac{0.15}{(0.5)(1)}=0.3.

Thus correlated outcomes produce correlated estimated means, even though each mean remains unbiased. Increasing the number of independent cases would reduce each covariance by the corresponding factor $1/n$ in this intercept-only model.

Worked Example: Predictor Geometry and Cross-Outcome Coefficients

Consider two predictors with design columns $x_1=(1,1,1,1)^\top$ and $x_2=(1,-1,0,0)^\top$. Suppose the error rows are independent bivariate normal vectors with $\Sigma=\begin{pmatrix}9&-3\\-3&4\end{pmatrix}$. Calculate the covariance and correlation between the two estimated outcome coefficients for each predictor.

The columns are orthogonal: $x_1^\top x_2=1-1+0+0=0$. Their squared lengths are $x_1^\top x_1=4$ and $x_2^\top x_2=2$. Therefore,

X^\top X=\begin{pmatrix}4&0\\0&2\end{pmatrix}, \qquad C=\begin{pmatrix}1/4&0\\0&1/2\end{pmatrix}.

For predictor 1, the two outcome-specific coefficient variances are $(1/4)(9)=2.25$ and $(1/4)(4)=1$, while their covariance is $(1/4)(-3)=-0.75$. Their correlation is $-0.75/(\sqrt{2.25}\sqrt{1})=-0.5$. For predictor 2, the variances are $(1/2)(9)=4.5$ and $(1/2)(4)=2$, and the covariance is $(1/2)(-3)=-1.5$; the correlation is $-1.5/(\sqrt{4.5}\sqrt{2})=-0.5$. The predictor estimates have zero covariance across predictors because $C_{12}=0$. Orthogonal predictor columns remove that source of coefficient covariance, but do not remove dependence between outcome-specific estimates for the same predictor.

What the Formula Says—and a Common Pitfall

For a fixed predictor $r$, if both outcome variances are positive, the correlation between the estimates for outcomes $j$ and $k$ is

\operatorname{Corr}(\hat B_{rj},\hat B_{rk}) =\frac{\Sigma_{jk}}{\sqrt{\Sigma_{jj}\Sigma_{kk}}}.

The design factor $C_{rr}$ cancels. Thus, for a given predictor, the correlation between its estimated coefficients across outcomes equals the within-case error correlation between those outcomes. Across different predictors, the design factor $C_{rs}$ also enters; it can make coefficient covariance positive, negative, or zero.

A common mistake is to reason that correlated outcomes make OLS coefficients biased. That conclusion does not follow: the expectation calculation uses $E[E]=0$, and correlation between outcomes does not change that condition. Correlation does affect the joint sampling variation, as the covariance formula shows. A separate mistake is to assume that zero covariance always means independence. The derivation establishes zero covariance when a factor in the product is zero; it does not establish independence without additional distributional assumptions.

These results also clarify the scope of the assumptions. If $X$ is random, the same calculations apply conditionally on $X$ when the conditional error mean is zero, the conditional row covariances are $\Sigma$, and the conditional rows are independent. Merely having zero unconditional error means is not enough when the errors may depend on the design. The full-column-rank condition is needed for the stated inverse formula.

Key takeaway: With a full-rank fixed design and independent, zero-mean error rows sharing covariance $\Sigma$, multivariate OLS is unbiased and its coefficient covariance is $\Sigma\otimes(X^\top X)^{-1}$. The design governs variation across predictors; the outcome covariance governs variation across responses.

Check Your Understanding

Use the model conditions and coefficient covariance formula to answer the following.

Which assumption gives $E[\hat B]=B$, and which assumptions are additionally used to derive its covariance?
If $C_{rr}=0.4$, $\Sigma_{jj}=5$, and $\Sigma_{jk}=-1$, find $\operatorname{Var}(\hat B_{rj})$ and $\operatorname{Cov}(\hat B_{rj},\hat B_{rk})$.
When does the covariance between two coefficients for different predictors vanish, and why does that not by itself imply independence?
For a fixed predictor, what is the correlation between its two outcome-specific coefficient estimates when the corresponding error correlation is $0.25$?
Explain why correlated outcomes can change coefficient uncertainty without making OLS biased.