The General Multivariate Linear Model

Why a Multivariate Linear Model?

In many studies, each case contributes several related outcomes: a patient may have multiple clinical measurements, or a school may report several performance measures. At the same time, the expected outcomes may depend on predictors such as treatment group, age, or baseline characteristics. The general multivariate linear model expresses these relationships together, while allowing each outcome to have its own set of regression coefficients.

The matrix notation is compact, but its dimensions carry the model’s meaning. The rows correspond to cases, the columns of the response matrix to outcomes, and the columns of the design matrix to predictor terms. The coefficient matrix connects those predictor terms to the outcomes; the error matrix records deviations from the modeled values.

Definition: The general multivariate linear model is

$$ Y=XB+U, $$

where \(Y\) is an \(n\times p\) response matrix, \(X\) is an \(n\times q\) design matrix, \(B\) is a \(q\times p\) coefficient matrix, and \(U\) is an \(n\times p\) error matrix. Thus \(XB\) is also \(n\times p\), and the equation is dimensionally defined.

Reading the Model Entry by Entry

Let \(y_{ij}\) denote the response for case \(i\) on outcome \(j\), let \(x_{ik}\) be the value of design term \(k\) for case \(i\), and let \(b_{kj}\) be the coefficient connecting design term \(k\) to outcome \(j\). The matrix equation says, for every case \(i\) and outcome \(j\),

y_{ij}=\sum_{k=1}^{q}x_{ik}b_{kj}+u_{ij}, \qquad i=1,\ldots,n,\quad j=1,\ldots,p.

The same \(X\) is used for all \(p\) outcomes in this formulation, while the columns of \(B\) allow the predictor effects to differ by outcome. A row of \(Y\) is the vector of observed outcomes for one case. A row of \(U\) is that case’s vector of errors. The columns of \(X\) are design terms, which may include an intercept, measured predictors, or indicators for groups.

Worked Example: A Shared Design With Two Outcomes

Suppose \(n=4\) patients have two outcomes, and the design includes an intercept and one predictor \(x\). Then \(p=2\) and \(q=2\), so \(Y\) and \(U\) are \(4\times2\), \(X\) is \(4\times2\), and \(B\) is \(2\times2\). For patient \(i\), the model gives

\begin{aligned} y_{i1}&=b_{11}+x_i b_{21}+u_{i1},\\ y_{i2}&=b_{12}+x_i b_{22}+u_{i2}. \end{aligned}

For example, if \(x_3=5\), the modeled values for patient 3 are \(b_{11}+5b_{21}\) and \(b_{12}+5b_{22}\). The coefficients \(b_{21}\) and \(b_{22}\) need not be equal: the predictor can have different effects on the two outcomes. The error terms \(u_{31}\) and \(u_{32}\) are the deviations for this patient on the respective outcomes.

Worked Example: Checking Matrix Dimensions

A study observes \(n=60\) people, records \(p=3\) responses per person, and uses \(q=5\) design terms, including an intercept. The required dimensions are

Y:60\times3,\qquad X:60\times5,\qquad B:5\times3,\qquad U:60\times3.

Multiplication gives \(XB\) dimensions \((60\times5)(5\times3)=60\times3\), matching \(Y\) and \(U\). Each column of \(B\) contains the five coefficients for one response. There are \(5\cdot3=15\) coefficient entries in total. This count alone does not guarantee they can all be separately identified: the columns of \(X\) must also contain enough independent information.

Mean Structure and Error Assumptions

The equation \(Y=XB+U\) defines a decomposition, but statistical interpretation requires assumptions about the errors. One common formulation treats \(X\) as fixed. Another treats it as random and interprets conditions given the observed design. In either case, a basic mean-model assumption is that errors have mean zero, conditional on the design.

Mean condition: When the expectations exist, the zero conditional mean assumption is \(E[U\mid X]=0\), interpreted as \(E[U]=0\) for a fixed design. It implies \(E[Y\mid X]=XB\). For estimation by squared error, finite second moments are also typically assumed so that the expected squared errors are finite.

This mean condition does not specify the covariance structure of the errors. In particular, it does not by itself assert that different cases are independent, that responses within a case are uncorrelated, or that all errors have the same variance. Such assumptions may be needed for particular standard errors or distributional conclusions, but they are not required simply to write the model or define least squares.

The coefficient matrix is a meaningful target only if the design distinguishes its entries. The relevant design condition is full column rank: the \(q\) columns of \(X\) are linearly independent, or equivalently \(\operatorname{rank}(X)=q\). Necessarily \(q\le n\). This is an identifiability condition, not a claim about the probability distribution of \(U\).

Worked Example: A Redundant Design Column

Suppose a design has an intercept column and two predictor columns, with the second predictor equal to twice the first for every case. If the first predictor column is \(x\) and the intercept column is \(\mathbf{1}\), then

X=\begin{bmatrix}\mathbf{1}&x&2x\end{bmatrix}.

The third column is twice the second, so the columns are dependent and \(\operatorname{rank}(X)<3\). For one response, the linear predictor is \(b_1\mathbf{1}+b_2x+b_3(2x)=b_1\mathbf{1}+(b_2+2b_3)x\). Only the combination \(b_2+2b_3\) affects the fitted mean. For instance, \((b_2,b_3)=(1,0)\) and \((b_2,b_3)=(3,-1)\) both give that combination equal to 1. The data design cannot distinguish those coefficient pairs.

Least-Squares Estimation and Rank

For observed \(Y\) and \(X\), least squares seeks a coefficient matrix whose fitted matrix \(XB\) is close to \(Y\), measured by the sum of squared entries of the residual matrix. Define the criterion

Q(B)=\sum_{i=1}^{n}\sum_{j=1}^{p}(y_{ij}-(XB)_{ij})^2 =\|Y-XB\|_F^2,

where \(\|A\|_F^2\) is the sum of squared entries of \(A\). A least-squares estimate is any matrix \(\hat B\) minimizing \(Q(B)\). No error independence or normality assumption is needed to define this optimization problem.

Theorem: For any finite matrices \(Y\) and \(X\) with compatible dimensions, a least-squares minimizer exists. It is unique if and only if \(\operatorname{rank}(X)=q\). If \(X\) is rank deficient, every minimizer has at least one distinct coefficient matrix with the same fitted matrix and the same criterion value.

Proof. To establish existence, consider one column \(y\) of \(Y\). The possible fitted columns \(Xb\) form the column space \(C\) of \(X\), a finite-dimensional subspace of Euclidean space and therefore a closed set. Let \(d=\inf_{z\in C}\|y-z\|\), and choose \(z_m\in C\) with \(\|y-z_m\|\to d\). The sequence is bounded because \(\|z_m\|\le \|y\|+\|y-z_m\|\), and the right-hand side is bounded. A bounded sequence in a finite-dimensional space has a convergent subsequence. Its limit belongs to \(C\), since \(C\) is closed, and continuity of the norm shows that this limit attains distance \(d\). Thus a best fitted column exists and has at least one coefficient vector \(b\) with \(Xb=z\). Applying this argument to each of the finitely many columns of \(Y\) gives a coefficient matrix minimizing \(Q(B)\).

If \(\operatorname{rank}(X)=q\), then \(XD\ne0\) for every nonzero \(q\times p\) matrix \(D\): at least one column of \(D\) is nonzero, and \(X\) has trivial null space. For any \(B_1\ne B_2\), set \(D=B_1-B_2\). The squared-error criterion satisfies

Q\!\left(\frac{B_1+B_2}{2}\right) =\frac{Q(B_1)+Q(B_2)}{2}-\frac{\|X(B_1-B_2)\|_F^2}{4}.

The final term is strictly positive when \(B_1\ne B_2\). Therefore two distinct matrices cannot both minimize \(Q\), so the minimizer is unique.

If \(\operatorname{rank}(X)<q\), there is a nonzero vector \(v\) such that \(Xv=0\). Let \(B\) be any minimizer, whose existence was just established, and choose any nonzero vector \(a\in\mathbb{R}^{p}\). Then \(va^T\ne0\), and

$$ X(B+va^T)=XB+(Xv)a^T=XB. $$

Hence \(B+va^T\) is distinct from \(B\), has the same fitted matrix, and gives the same value of \(Q\). The minimizer is not unique. This proves both directions. \(\square\)

What the Conditions Do—and Do Not—Guarantee

Full column rank makes the coefficient matrix identifiable within the stated linear mean model and gives a unique least-squares solution. It does not guarantee that the linear mean model is scientifically appropriate, that predictions are accurate, or that a particular inferential procedure has its claimed distribution. Those conclusions require additional assumptions suited to the goal of the analysis.

The zero conditional mean condition gives the model its mean interpretation: the conditional mean response is \(XB\). If it fails, the fitted linear structure need not represent the conditional mean, even when least squares can still be computed. Conversely, dependent or correlated errors do not make the matrix equation invalid; they affect the analysis of uncertainty and may require methods beyond ordinary least squares.

Key takeaway: The model’s dimensions encode how cases, outcomes, design terms, and coefficients align. A zero conditional error mean supports the mean interpretation, while full column rank of \(X\) is the condition that allows the full coefficient matrix to be uniquely identified by least squares.

Common Pitfall

A common mistake is to conclude that having at least as many cases as design terms, \(n\ge q\), guarantees unique coefficient estimates. That is not sufficient: two or more design columns may still be linearly dependent. The correct condition is \(\operatorname{rank}(X)=q\), which requires both \(n\ge q\) and independence of all \(q\) columns.

A second mistake is to treat the zero-mean condition as if it implied independent errors or uncorrelated outcomes. It does not. The condition \(E[U\mid X]=0\) concerns the conditional mean only; covariance and independence assumptions must be stated separately when needed.

Check Your Understanding

Use the matrix model and conditions from this tutorial to answer the following.

If there are 40 cases, 4 outcomes, and 6 design terms, what are the dimensions of \(Y\), \(X\), \(B\), and \(U\)?
In the entrywise model, what does \(b_{kj}\) describe?
Does \(E[U\mid X]=0\) imply that errors for different outcomes are uncorrelated? Explain.
Why is \(n\ge q\) alone insufficient to guarantee unique least-squares coefficients?
If \(Xv=0\) for a nonzero vector \(v\), explain how a least-squares minimizer can be changed without changing its fitted matrix.