Multivariate Linear Models: Zero to Hero Tutorial

This comprehensive tutorial takes you from the foundational concepts of Multivariate Linear Models (MLM) all the way through advanced estimation, hypothesis testing, model diagnostics, coefficient interpretation, and practical usage within the DataStatPro application. Whether you are a complete beginner or an experienced analyst, this guide is structured to build your understanding step by step.

Table of Contents

1. Prerequisites and Background Concepts
2. What are Multivariate Linear Models?
3. The Mathematical Framework
4. Model Specification and Design Matrices
5. Assumptions of Multivariate Linear Models
6. Parameter Estimation
7. Hypothesis Testing and Inference
8. Multivariate Test Statistics
9. Effect Size Measures
10. Model Fit and Evaluation
11. Model Diagnostics and Residuals
12. Interpretation of Coefficients
13. Variable Selection and Model Comparison
14. Special Cases and Connections to Other Methods
15. Using the Multivariate Linear Models Component
16. Computational and Formula Details
17. Worked Examples
18. Common Mistakes and How to Avoid Them
19. Troubleshooting
20. Quick Reference Cheat Sheet

1. Prerequisites and Background Concepts

Before diving into Multivariate Linear Models, it is helpful to be familiar with the following foundational concepts. Do not worry if you are not — each concept is briefly explained here.

1.1 Ordinary Least Squares (OLS) Regression

Ordinary Least Squares (OLS) regression models a single continuous response $y$ as a linear function of $p$ predictors $X_1, X_2, \dots, X_p$:

$y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \dots + \beta_p X_{ip} + \epsilon_i$

In matrix form:

$\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}, \quad \boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I})$

The OLS estimator minimises the sum of squared residuals:

$\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}$

Multivariate Linear Models extend this framework to the case where there are multiple response variables measured simultaneously.

1.2 Matrices and Linear Algebra

Key matrix operations used throughout this tutorial:

• Matrix transpose $\mathbf{A}^T$: Swaps rows and columns.
• Matrix inverse $\mathbf{A}^{-1}$: Exists when $\mathbf{A}$ is square and full rank; $\mathbf{A}\mathbf{A}^{-1} = \mathbf{I}$.
• Trace $\text{tr}(\mathbf{A}) = \sum_i a_{ii}$: Sum of diagonal elements.
• Determinant $|\mathbf{A}|$: Scalar summarising the matrix; zero if $\mathbf{A}$ is singular.
• Kronecker product $\mathbf{A} \otimes \mathbf{B}$: Block matrix of all pairwise products of elements of $\mathbf{A}$ with $\mathbf{B}$.
• Vec operator $\text{vec}(\mathbf{A})$: Stacks the columns of $\mathbf{A}$ into a single column vector.

1.3 The Multivariate Normal Distribution

The multivariate normal distribution $\mathcal{N}_q(\boldsymbol{\mu}, \boldsymbol{\Sigma})$ has density:

$f(\mathbf{y}) = \frac{1}{(2\pi)^{q/2}|\boldsymbol{\Sigma}|^{1/2}} \exp\left\{-\frac{1}{2}(\mathbf{y}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\mathbf{y}-\boldsymbol{\mu})\right\}$

Where $\boldsymbol{\mu}$ is the mean vector and $\boldsymbol{\Sigma}$ is the covariance matrix. MLMs assume that the response vectors follow this distribution conditionally on the predictors.

1.4 Covariance and Correlation Matrices

The covariance matrix $\boldsymbol{\Sigma}$ ($q \times q$) for $q$ response variables contains:

• Variances $\sigma_{jj} = \text{Var}(Y_j)$ on the diagonal.
• Covariances $\sigma_{jk} = \text{Cov}(Y_j, Y_k)$ off-diagonal.

The correlation matrix $\mathbf{R}$ standardises covariances:

$R_{jk} = \frac{\sigma_{jk}}{\sqrt{\sigma_{jj}\sigma_{kk}}}$

1.5 Eigenvalues and Eigenvectors

For a square matrix $\mathbf{A}$, eigenvalue-eigenvector pairs $(\lambda_s, \mathbf{v}_s)$ satisfy $\mathbf{A}\mathbf{v}_s = \lambda_s\mathbf{v}_s$. They are critical for computing multivariate test statistics and understanding the structure of relationships in MLMs.

1.6 The Hat Matrix

In OLS regression, the hat (projection) matrix:

$\mathbf{H} = \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T$

Projects the response vector $\mathbf{y}$ onto the column space of $\mathbf{X}$: $\hat{\mathbf{y}} = \mathbf{H}\mathbf{y}$. The hat matrix diagonal elements $h_{ii}$ measure the leverage of each observation. The same matrix plays the same role in MLMs, now projecting each response variable.

2. What are Multivariate Linear Models?

2.1 The Core Idea

A Multivariate Linear Model (MLM) — also called multivariate multiple regression or the general multivariate linear model — simultaneously models the linear relationship between a set of predictor variables (independent variables) and multiple continuous response variables (dependent variables).

Where univariate multiple regression has a single response $\mathbf{y}$ ($n \times 1$), the multivariate linear model has a response matrix $\mathbf{Y}$ ($n \times q$) with $q \geq 2$ response variables measured on the same $n$ observations.

The model is:

$\mathbf{Y} = \mathbf{X}\mathbf{B} + \mathbf{E}$

Where:

• $\mathbf{Y}$ ($n \times q$) = matrix of $q$ response variables.
• $\mathbf{X}$ ($n \times (p+1)$) = design matrix of predictors (including a column of ones for the intercept).
• $\mathbf{B}$ ($(p+1) \times q$) = matrix of regression coefficients (one row per predictor, one column per response).
• $\mathbf{E}$ ($n \times q$) = matrix of errors.

2.2 Why Use MLM Instead of Separate Regressions?

A natural question is: "Why not simply run $q$ separate univariate regressions, one for each response?" There are several compelling reasons to prefer the multivariate approach:

2.3 The Generality of the Multivariate Linear Model

The multivariate linear model is remarkably general. Many seemingly distinct statistical procedures are special cases:

2.4 Real-World Applications

Multivariate Linear Models are used across a wide range of applied fields:

• Neuroscience: Predicting multiple brain region activation levels simultaneously from experimental conditions, participant characteristics, or cognitive task measures.
• Environmental Science: Modelling multiple water quality indicators (pH, dissolved oxygen, turbidity, nitrate) jointly as functions of environmental predictors (rainfall, temperature, land use).
• Finance: Predicting multiple asset returns simultaneously from macroeconomic factors (interest rates, inflation, GDP growth).
• Clinical Medicine: Relating treatment variables to multiple physiological endpoints (blood pressure, cholesterol, glucose, BMI) simultaneously.
• Educational Research: Modelling multiple academic outcomes (reading, mathematics, science, writing) as joint functions of student and school characteristics.
• Agricultural Science: Predicting multiple crop yield components (grain weight, harvest index, protein content) from soil, climate, and management variables.
• Marketing Research: Modelling multiple consumer attitude dimensions simultaneously from demographic and psychographic predictors.
• Psychometrics: Relating latent constructs to multiple observed scale scores while accounting for the shared measurement structure.

3. The Mathematical Framework

3.1 The Multivariate Linear Model

The core model is:

$\mathbf{Y}_{n \times q} = \mathbf{X}_{n \times (p+1)} \mathbf{B}_{(p+1) \times q} + \mathbf{E}_{n \times q}$

Each row $\mathbf{y}_i^T$ (observation $i$) satisfies:

$\mathbf{y}_i = \mathbf{B}^T \mathbf{x}_i + \boldsymbol{\epsilon}_i, \quad \boldsymbol{\epsilon}_i \sim \mathcal{N}_q(\mathbf{0}, \boldsymbol{\Sigma})$

Where:

• $\mathbf{x}_i$ ($p+1 \times 1$) = predictor vector for observation $i$ (including intercept).
• $\boldsymbol{\Sigma}$ ($q \times q$) = error covariance matrix (assumed common to all observations).
• Errors are independent across observations: $\text{Cov}(\boldsymbol{\epsilon}_i, \boldsymbol{\epsilon}_j) = \mathbf{0}$ for $i \neq j$.

Written differently, the distribution of each row is:

$\mathbf{y}_i^T \mid \mathbf{x}_i \sim \mathcal{N}_q\left(\mathbf{x}_i^T \mathbf{B},\ \boldsymbol{\Sigma}\right)$

3.2 The Coefficient Matrix $\mathbf{B}$

The coefficient matrix $\mathbf{B}$ has dimensions $(p+1) \times q$:

$\mathbf{B} = \begin{pmatrix} \beta_{01} & \beta_{02} & \cdots & \beta_{0q} \\ \beta_{11} & \beta_{12} & \cdots & \beta_{1q} \\ \vdots & \vdots & \ddots & \vdots \\ \beta_{p1} & \beta_{p2} & \cdots & \beta_{pq} \end{pmatrix}$

• The $j$-th column $\boldsymbol{\beta}_j = (\beta_{0j}, \beta_{1j}, \dots, \beta_{pj})^T$ contains the coefficients for the $j$-th response variable — identical to what you would get from running univariate OLS regression of $Y_j$ on $\mathbf{X}$.
• The $k$-th row (for $k \geq 1$) contains the effects of predictor $X_k$ on all $q$ response variables simultaneously — the multivariate coefficient vector for predictor $k$.

3.3 The Error Structure

The error matrix $\mathbf{E}$ has the following properties:

$E[\mathbf{E}] = \mathbf{0}_{n \times q}$

$E[\mathbf{e}_i \mathbf{e}_i^T] = \boldsymbol{\Sigma}, \quad E[\mathbf{e}_i \mathbf{e}_j^T] = \mathbf{0} \text{ for } i \neq j$

This means the errors are identically and independently distributed (i.i.d.) multivariate normal. Equivalently, using the vec operator:

$\text{vec}(\mathbf{E}) \sim \mathcal{N}_{nq}\left(\mathbf{0},\ \boldsymbol{\Sigma} \otimes \mathbf{I}_n\right)$

The Kronecker product $\boldsymbol{\Sigma} \otimes \mathbf{I}_n$ captures the covariance structure: observations are independent (block-diagonal structure from $\mathbf{I}_n$) but responses within each observation may be correlated (captured by $\boldsymbol{\Sigma}$).

3.4 The Conditional Mean and Variance

For a new observation with predictor vector $\mathbf{x}_{new}$:

$E[\mathbf{y}_{new} \mid \mathbf{x}_{new}] = \mathbf{B}^T \mathbf{x}_{new}$

$\text{Var}(\mathbf{y}_{new} \mid \mathbf{x}_{new}) = \boldsymbol{\Sigma}$

The predicted response vector is:

$\hat{\mathbf{y}}_{new} = \hat{\mathbf{B}}^T \mathbf{x}_{new}$

The predicted response matrix for all $n$ observations:

$\hat{\mathbf{Y}} = \mathbf{X}\hat{\mathbf{B}} = \mathbf{H}\mathbf{Y}$

Where $\mathbf{H} = \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T$ is the hat matrix.

3.5 The General Linear Hypothesis

The power of the multivariate linear model framework lies in the ability to test very general hypotheses of the form:

$H_0: \mathbf{C}\mathbf{B}\mathbf{M} = \mathbf{\Gamma}_0$

Where:

• $\mathbf{C}$ ($c \times (p+1)$) = hypothesis (contrast) matrix for predictors — specifies which linear combinations of rows of $\mathbf{B}$ are being tested.
• $\mathbf{M}$ ($q \times m$) = response transformation matrix — specifies which linear combinations of columns of $\mathbf{B}$ (i.e., response variables or contrasts among responses) are being tested.
• $\mathbf{\Gamma}_0$ ($c \times m$) = hypothesised value matrix (usually $\mathbf{0}$ for testing no effect).

This general framework subsumes all the special cases mentioned in Section 2.3.

Examples of $\mathbf{C}\mathbf{B}\mathbf{M} = \mathbf{0}$ hypotheses:

4. Model Specification and Design Matrices

4.1 The Design Matrix $\mathbf{X}$

The design matrix $\mathbf{X}$ ($n \times (p+1)$) encodes all predictor information. The first column is typically a column of ones for the intercept:

$\mathbf{X} = \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1p} \\ 1 & x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{np} \end{pmatrix}$

The design matrix can accommodate:

• Continuous predictors: Entered directly as columns.
• Categorical predictors: Encoded as dummy variables ($k-1$ columns for $k$ categories).
• Interaction terms: Products of two or more predictor columns.
• Polynomial terms: Powers of continuous predictors ($X^2$, $X^3$, etc.).
• Mixed designs: Combinations of continuous and categorical predictors.

4.2 Dummy Coding for Categorical Predictors

For a categorical predictor with $k$ categories, create $k-1$ dummy variables using a reference category. For a three-level factor (A, B, C) with A as reference:

The coefficient for $D_B$ in column $j$ of $\hat{\mathbf{B}}$ represents the difference in the mean of $Y_j$ between groups B and A, holding other predictors constant.

4.3 The Response Matrix $\mathbf{Y}$

The response matrix $\mathbf{Y}$ ($n \times q$) contains all $q$ response variables:

$\mathbf{Y} = \begin{pmatrix} y_{11} & y_{12} & \cdots & y_{1q} \\ y_{21} & y_{22} & \cdots & y_{2q} \\ \vdots & \vdots & \ddots & \vdots \\ y_{n1} & y_{n2} & \cdots & y_{nq} \end{pmatrix}$

Each column is a separate response variable. Each row is the complete response profile for one observation.

4.4 Centering and Scaling

Centering predictors (subtracting their means) is generally recommended to:

• Make the intercept interpretable (the predicted response at the mean of all predictors).
• Reduce multicollinearity between main effects and interactions.
• Improve numerical stability of matrix computations.

Standardising predictors (centering + scaling to unit variance) additionally makes regression coefficients comparable across predictors with different units and scales — the standardised coefficients represent the change in the response (in standard deviation units of the response) for a one-standard-deviation change in each predictor.

4.5 Interaction Terms

An interaction term $X_j \times X_k$ is formed as the element-wise product of the two predictor columns:

$\mathbf{x}_{jk} = \mathbf{x}_j \odot \mathbf{x}_k$

Including an interaction in the model allows the effect of $X_j$ on the response vector to depend on the value of $X_k$ (and vice versa). For categorical × continuous interactions, the interaction allows the regression slope to differ across groups.

5. Assumptions of Multivariate Linear Models

5.1 Linearity

Assumption: The relationship between each predictor $X_k$ and each response $Y_j$ is linear (on the scale of the model), holding all other predictors constant:

$E[Y_{ij} \mid \mathbf{x}_i] = \beta_{0j} + \beta_{1j}X_{i1} + \dots + \beta_{pj}X_{ip}$

How to check: Partial residual plots (component-plus-residual plots) for each predictor-response combination; scatter plots of residuals against each predictor.

Consequences of violation: Biased and inconsistent coefficient estimates; poor prediction.

Remedies: Polynomial terms, log transformations of predictors or responses, spline terms, generalised additive models.

5.2 Multivariate Normality of Errors

Assumption: The error vectors $\boldsymbol{\epsilon}_i$ follow a multivariate normal distribution:

$\boldsymbol{\epsilon}_i \sim \mathcal{N}_q(\mathbf{0}, \boldsymbol{\Sigma})$

How to check:

• Mardia's tests of multivariate skewness and kurtosis on residuals.
• Q-Q plots of squared Mahalanobis distances of residuals vs. $\chi^2_q$ quantiles.
• Univariate normality checks (histograms, Q-Q plots, Shapiro-Wilk) on each column of $\hat{\mathbf{E}}$.

Consequences of violation: OLS estimates remain unbiased and consistent (Gauss-Markov theorem extends to MLM) but inference (hypothesis tests, confidence intervals) based on the multivariate normal assumption is affected, particularly with small samples.

Remedies: Transformations (log, Box-Cox) of skewed responses; robust estimation; bootstrap inference.

5.3 Independence of Observations

Assumption: The error vectors $\boldsymbol{\epsilon}_i$ are independent across observations:

$\text{Cov}(\boldsymbol{\epsilon}_i, \boldsymbol{\epsilon}_j) = \mathbf{0} \quad \text{for } i \neq j$

How to check: Consider the study design. Look for temporal autocorrelation (Durbin-Watson statistics per response), spatial correlation, or clustering structure.

Consequences of violation: Biased standard errors, invalid hypothesis tests, false significance.

Remedies: Mixed models (multilevel MLM) for clustered data; GEE for repeated measures; time-series corrections for temporal dependence.

5.4 Homoscedasticity (Constant Error Covariance)

Assumption: The error covariance matrix $\boldsymbol{\Sigma}$ is the same for all observations — it does not depend on the values of the predictors or the fitted values:

$\text{Var}(\boldsymbol{\epsilon}_i) = \boldsymbol{\Sigma} \quad \text{for all } i$

How to check:

• Plots of squared residuals vs. fitted values for each response.
• Levene's test or Bartlett's test for each response variable.
• Box's M test if comparing across groups defined by a categorical predictor.

Consequences of violation: OLS estimates remain unbiased but are no longer BLUE (Best Linear Unbiased Estimator); standard errors are biased.

Remedies: Weighted least squares (if the heteroscedasticity pattern is known); heteroscedasticity-consistent (HC) standard errors; transformations of responses.

5.5 No Perfect Multicollinearity Among Predictors

Assumption: No predictor variable is a perfect linear combination of other predictor variables. Equivalently, $\mathbf{X}^T\mathbf{X}$ must be invertible (full rank).

How to check: Variance Inflation Factor (VIF) for each predictor: $VIF_k = 1/(1 - R^2_k)$ where $R^2_k$ is the R² from regressing $X_k$ on all other predictors. VIF > 10 is a common threshold for concern.

Consequences of violation: $(\mathbf{X}^T\mathbf{X})^{-1}$ does not exist; coefficient estimates are undefined or numerically unstable with inflated standard errors.

Remedies: Remove redundant predictors; use ridge regression or PLS; combine highly correlated predictors into a composite.

5.6 No Influential Outliers

Assumption: No single observation has undue influence on the estimated coefficient matrix $\hat{\mathbf{B}}$.

How to check: Cook's distance (multivariate extension), hat matrix diagonal $h_{ii}$, DFFITS, standardised residuals.

Consequences of violation: Estimated coefficients may be heavily distorted by a single atypical observation.

Remedies: Investigate outlying observations for data errors; use robust regression; report sensitivity analyses with and without influential points.

5.7 Sufficient Sample Size

Recommendation: For MLM to be reliable:

• $n > p + q + 1$ (absolute minimum to estimate $\hat{\mathbf{B}}$ and $\hat{\boldsymbol{\Sigma}}$).
• $n \geq 10(p + 1)$ (rule of thumb for adequate precision of coefficient estimates).
• The within-group covariance matrix $\hat{\boldsymbol{\Sigma}}$ requires $n - p - 1 > q$ for invertibility.

6. Parameter Estimation

6.1 The Ordinary Least Squares (OLS) Estimator

The OLS estimator for $\mathbf{B}$ minimises the total sum of squared residuals simultaneously across all response variables:

$\hat{\mathbf{B}} = \arg\min_{\mathbf{B}} \text{tr}\left[(\mathbf{Y} - \mathbf{X}\mathbf{B})^T(\mathbf{Y} - \mathbf{X}\mathbf{B})\right]$

Taking the matrix derivative and setting to zero:

$\frac{\partial}{\partial \mathbf{B}}\text{tr}\left[(\mathbf{Y} - \mathbf{X}\mathbf{B})^T(\mathbf{Y} - \mathbf{X}\mathbf{B})\right] = -2\mathbf{X}^T(\mathbf{Y} - \mathbf{X}\mathbf{B}) = \mathbf{0}$

Solving the normal equations $\mathbf{X}^T\mathbf{X}\hat{\mathbf{B}} = \mathbf{X}^T\mathbf{Y}$:

$\hat{\mathbf{B}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{Y}$

Critically: This is simply the matrix of coefficients obtained by running $q$ separate OLS regressions — one for each column of $\mathbf{Y}$. The multivariate OLS estimator is column-by-column identical to the univariate OLS estimator applied $q$ times.

6.2 Properties of the OLS Estimator

Under the standard assumptions:

Unbiasedness: $E[\hat{\mathbf{B}}] = \mathbf{B}$

Covariance structure:

$\text{Cov}(\hat{\boldsymbol{\beta}}_j, \hat{\boldsymbol{\beta}}_k) = \sigma_{jk}(\mathbf{X}^T\mathbf{X})^{-1}$

Where $\sigma_{jk}$ is the $(j,k)$ element of $\boldsymbol{\Sigma}$ (the error covariance between responses $j$ and $k$). In full:

$\text{Cov}(\text{vec}(\hat{\mathbf{B}})) = \boldsymbol{\Sigma} \otimes (\mathbf{X}^T\mathbf{X})^{-1}$

Gauss-Markov (multivariate): Among all linear unbiased estimators, $\hat{\mathbf{B}}$ has the minimum covariance matrix (in the Loewner partial order) — it is BLUE.

Maximum Likelihood Equivalence: Under multivariate normality, the OLS estimator is also the maximum likelihood estimator (MLE) of $\mathbf{B}$.

6.3 Estimation of the Error Covariance Matrix

The unbiased estimator of the error covariance matrix $\boldsymbol{\Sigma}$ is:

$\hat{\boldsymbol{\Sigma}} = \frac{\hat{\mathbf{E}}^T\hat{\mathbf{E}}}{n - p - 1} = \frac{(\mathbf{Y} - \hat{\mathbf{Y}})^T(\mathbf{Y} - \hat{\mathbf{Y}})}{n - p - 1}$

Where $\hat{\mathbf{E}} = \mathbf{Y} - \hat{\mathbf{Y}} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$ is the matrix of OLS residuals.

The MLE of $\boldsymbol{\Sigma}$ (biased but with maximum likelihood properties) is:

$\tilde{\boldsymbol{\Sigma}}_{MLE} = \frac{\hat{\mathbf{E}}^T\hat{\mathbf{E}}}{n}$

The unbiased estimator $\hat{\boldsymbol{\Sigma}}$ is generally preferred for inference.

Degrees of freedom: The residual SSCP matrix $\hat{\mathbf{E}}^T\hat{\mathbf{E}}$ has $n - p - 1$ degrees of freedom. Invertibility requires $n - p - 1 > q$.

6.4 Seemingly Unrelated Regression (SUR) and GLS

When the predictor matrices differ across response variables (i.e., each response has its own set of predictors), standard OLS on each response separately is not the most efficient estimator. Seemingly Unrelated Regression (SUR) applies Generalised Least Squares (GLS) using the full covariance structure:

$\hat{\mathbf{B}}_{GLS} = \left[\mathbf{X}^T(\hat{\boldsymbol{\Sigma}}^{-1} \otimes \mathbf{I}_n)\mathbf{X}\right]^{-1}\mathbf{X}^T(\hat{\boldsymbol{\Sigma}}^{-1} \otimes \mathbf{I}_n)\text{vec}(\mathbf{Y})$

When all responses share the same predictor matrix $\mathbf{X}$ (the standard MLM case), SUR and OLS give identical estimates. The GLS advantage only manifests when predictor matrices differ across equations.

6.5 The Maximum Likelihood Estimator

The log-likelihood for the multivariate linear model is:

$\ell(\mathbf{B}, \boldsymbol{\Sigma}) = -\frac{n}{2}\ln|\boldsymbol{\Sigma}| - \frac{1}{2}\text{tr}\left[\boldsymbol{\Sigma}^{-1}(\mathbf{Y} - \mathbf{X}\mathbf{B})^T(\mathbf{Y} - \mathbf{X}\mathbf{B})\right] + \text{const}$

Maximising over $\mathbf{B}$ gives $\hat{\mathbf{B}}_{MLE} = \hat{\mathbf{B}}_{OLS}$ (as stated above). Maximising over $\boldsymbol{\Sigma}$ gives $\tilde{\boldsymbol{\Sigma}}_{MLE} = \hat{\mathbf{E}}^T\hat{\mathbf{E}}/n$.

The maximised log-likelihood is:

$\ell(\hat{\mathbf{B}}, \tilde{\boldsymbol{\Sigma}}_{MLE}) = -\frac{n}{2}\left[q\ln(2\pi) + \ln|\tilde{\boldsymbol{\Sigma}}_{MLE}| + q\right]$

7. Hypothesis Testing and Inference

7.1 The General Linear Hypothesis Framework

The unified framework for hypothesis testing in MLM is the General Linear Hypothesis:

$H_0: \mathbf{C}\mathbf{B}\mathbf{M} = \mathbf{\Gamma}_0$

For testing $H_0: \mathbf{C}\mathbf{B}\mathbf{M} = \mathbf{0}$ (the most common case), the Hypothesis SSCP matrix is:

$\mathbf{H} = \mathbf{M}^T\hat{\mathbf{B}}^T\mathbf{C}^T\left[\mathbf{C}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{C}^T\right]^{-1}\mathbf{C}\hat{\mathbf{B}}\mathbf{M}$

The Error SSCP matrix (the same for all hypotheses):

$\mathbf{E} = \mathbf{M}^T\hat{\mathbf{E}}^T\hat{\mathbf{E}}\mathbf{M} = \mathbf{M}^T(\mathbf{Y} - \hat{\mathbf{Y}})^T(\mathbf{Y} - \hat{\mathbf{Y}})\mathbf{M}$

When $\mathbf{M} = \mathbf{I}_q$, $\mathbf{E} = \hat{\mathbf{E}}^T\hat{\mathbf{E}}$.

Test statistics are based on the eigenvalues of $\mathbf{E}^{-1}\mathbf{H}$ (see Section 8).

7.2 Testing Individual Predictor Effects (Row of $\mathbf{B}$)

To test whether predictor $X_k$ has any effect on the set of responses ($H_0: \boldsymbol{\beta}_{k\cdot} = \mathbf{0}^T$, the $k$-th row of $\mathbf{B}$):

Set $\mathbf{C} = \mathbf{e}_k^T$ (row vector with 1 in position $k$, zeros elsewhere) and $\mathbf{M} = \mathbf{I}_q$.

This produces a multivariate test of whether $X_k$ simultaneously predicts zero for all $q$ responses. The four multivariate test statistics (Section 8) are applied.

7.3 Testing a Subset of Predictors

To test whether a subset of $c$ predictors jointly contribute to the model ($H_0: \mathbf{C}\mathbf{B} = \mathbf{0}$ for a $c \times (p+1)$ matrix $\mathbf{C}$):

The hypothesis SSCP matrix with $df_H = c$ is compared to the error SSCP matrix with $df_E = n - p - 1$.

7.4 Testing Specific Linear Combinations of Responses

To test whether a predictor has differential effects on different responses (e.g., $H_0: \beta_{k1} = \beta_{k2}$):

Set $\mathbf{C} = \mathbf{e}_k^T$ and $\mathbf{M} = (1, -1)^T$.

The transformed response is $Y_1 - Y_2$; the test becomes a univariate $t$-test on the difference.

7.5 Univariate Tests Within the MLM Framework

For each response variable $Y_j$ individually, the standard univariate $F$-test for predictor $X_k$ is obtained by setting $\mathbf{M} = \mathbf{e}_j$ and applying the Wald test:

$F_{kj} = \frac{\hat{\beta}_{kj}^2 / \hat{\sigma}_{jj}^{(kk)}}{1} \sim F_{1, n-p-1}$

Where $\hat{\sigma}_{jj}^{(kk)} = \hat{\sigma}_{jj} [(\mathbf{X}^T\mathbf{X})^{-1}]_{kk}$ is the estimated variance of $\hat{\beta}_{kj}$.

7.6 Wald Tests for Individual Coefficients

For a single coefficient $\beta_{kj}$ (predictor $k$, response $j$):

$t_{kj} = \frac{\hat{\beta}_{kj}}{SE(\hat{\beta}_{kj})} \sim t_{n-p-1}$

Where:

$SE(\hat{\beta}_{kj}) = \sqrt{\hat{\sigma}_{jj} \cdot \left[(\mathbf{X}^T\mathbf{X})^{-1}\right]_{kk}}$

A $(1-\alpha) \times 100\%$ confidence interval for $\beta_{kj}$:

$\hat{\beta}_{kj} \pm t_{\alpha/2,\, n-p-1} \times SE(\hat{\beta}_{kj})$

7.7 Simultaneous Confidence Regions

A key advantage of MLM is the ability to form joint confidence regions for multivariate coefficient vectors. The $100(1-\alpha)\%$ joint confidence region for the $k$-th row $\boldsymbol{\beta}_{k\cdot}$ of $\mathbf{B}$ (all $q$ responses simultaneously) is an ellipsoid:

$\left(\hat{\boldsymbol{\beta}}_{k\cdot} - \boldsymbol{\beta}_{k\cdot}\right)^T \left[\hat{\boldsymbol{\Sigma}} \cdot \left(\mathbf{X}^T\mathbf{X}\right)^{-1}_{kk}\right]^{-1} \left(\hat{\boldsymbol{\beta}}_{k\cdot} - \boldsymbol{\beta}_{k\cdot}\right) \leq q F_{\alpha,\, q,\, n-p-q}$

7.8 Testing the Overall Model

The overall model tests whether any predictor (excluding the intercept) has a significant effect on any response:

$H_0: \mathbf{B}_{-0} = \mathbf{0}_{p \times q}$

Where $\mathbf{B}_{-0}$ is $\mathbf{B}$ with the intercept row removed. This is tested with $\mathbf{C} = [\mathbf{0}_p \mid \mathbf{I}_p]$ (omitting the intercept column) and $\mathbf{M} = \mathbf{I}_q$.

8. Multivariate Test Statistics

All multivariate tests in MLM are based on the eigenvalues $\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_{s^*}$ of $\mathbf{E}^{-1}\mathbf{H}$, where $s^* = \min(c, m)$, $c$ = rank of $\mathbf{C}$ (number of hypothesis df), $m$ = number of transformed responses (columns of $\mathbf{M}$).

8.1 Wilks' Lambda ($\Lambda^*$)

$\Lambda^* = \frac{|\mathbf{E}|}{|\mathbf{H} + \mathbf{E}|} = \prod_{s=1}^{s^*} \frac{1}{1 + \lambda_s}$

Range: $0 \leq \Lambda^* \leq 1$; smaller values indicate stronger effects.

F-approximation (Rao):

$F = \frac{1 - \Lambda^{*1/t}}{\Lambda^{*1/t}} \cdot \frac{df_2}{df_1}$

Where:

$t = \sqrt{\frac{m^2 c^2 - 4}{m^2 + c^2 - 5}}, \quad df_1 = mc, \quad df_2 = t\left(df_E - \frac{m - c + 1}{2}\right) - \frac{mc - 2}{2}$

Exact F when $m = 1$, $m = 2$, $c = 1$, or $c = 2$.

8.2 Pillai's Trace ($V$)

$V = \text{tr}\left[\mathbf{H}(\mathbf{H}+\mathbf{E})^{-1}\right] = \sum_{s=1}^{s^*}\frac{\lambda_s}{1+\lambda_s}$

Range: $0 \leq V \leq s^*$.

F-approximation:

$F = \frac{V/s^*}{(s^* - V)/s^*} \cdot \frac{2n_Y + s^* + 1}{2n_X + s^* + 1} \sim F_{s^*(2n_X+s^*+1),\, s^*(2n_Y+s^*+1)}$

Where $n_X = (df_E - m - 1)/2$ and $n_Y = (c - m - 1)/2$.

Most robust to violations of assumptions.

8.3 Hotelling-Lawley Trace ($U$)

$U = \text{tr}(\mathbf{E}^{-1}\mathbf{H}) = \sum_{s=1}^{s^*}\lambda_s$

Range: $0 \leq U < \infty$.

F-approximation:

$F = \frac{U \cdot df_E(df_E - m - 1)}{s^* \cdot m \cdot (df_E + c)(df_E - m - 1) - 2U} \quad \text{(approximate)}$

Most powerful when effects are concentrated on a single dimension.

8.4 Roy's Largest Root ($\theta$)

$\theta = \frac{\lambda_1}{1 + \lambda_1}$

Range: $0 \leq \theta \leq 1$.

Roy's provides an upper bound to the $p$-value distribution. Exact tables are needed for precise $p$-values.

8.5 Comparison of the Four Statistics in MLM Context

8.6 Likelihood Ratio Test

An alternative formulation: Under $H_0: \mathbf{C}\mathbf{B}\mathbf{M} = \mathbf{0}$, the likelihood ratio statistic is:

$\Lambda_{LR} = \left(\frac{|\mathbf{E}|}{|\mathbf{H}+\mathbf{E}|}\right)^{n/2} = (\Lambda^*)^{n/2}$

The corresponding test statistic:

$\chi^2 \approx -\left[n - p - 1 - \frac{m + c + 1}{2}\right]\ln(\Lambda^*)$

Asymptotically $\sim \chi^2_{mc}$ under $H_0$. This approximation is Bartlett's correction to the likelihood ratio test.

9. Effect Size Measures

9.1 Multivariate Eta-Squared ($\eta^2_p$)

The most widely reported multivariate effect size for each hypothesis test:

From Wilks' Lambda:

$\eta^2_p = 1 - \Lambda^{*1/t}$

Where $t$ is defined as in Section 8.1.

From Pillai's Trace:

$\eta^2_p \approx \frac{V}{s^*}$

Benchmarks:

9.2 Multivariate Omega-Squared ($\omega^2$)

A less biased (corrected) multivariate effect size:

$\omega^2 = \frac{df_H(F - 1)}{df_H(F - 1) + n}$

Applied to the $F$-approximation from the multivariate tests. Can be computed separately from each test statistic's $F$ approximation.

9.3 Univariate Effect Sizes per Response

For each response variable $Y_j$, report the univariate $R^2$:

$R^2_j = 1 - \frac{SS_{res,j}}{SS_{tot,j}} = 1 - \frac{\hat{\mathbf{e}}_j^T\hat{\mathbf{e}}_j}{\sum_i(y_{ij}-\bar{y}_j)^2}$

And univariate partial $\eta^2_p$ for each predictor-response combination:

$\eta^2_{p,kj} = \frac{SS_{k,j}}{SS_{k,j} + SS_{res,j}}$

9.4 Standardised Coefficients

Standardised regression coefficients ($\beta^*_{kj}$) are obtained by standardising both predictors and responses before fitting:

$\beta^*_{kj} = \hat{\beta}_{kj} \cdot \frac{s_{X_k}}{s_{Y_j}}$

Where $s_{X_k}$ and $s_{Y_j}$ are the standard deviations of predictor $k$ and response $j$. Standardised coefficients represent the change in $Y_j$ (in SD units) for a one-SD change in $X_k$, facilitating comparison across predictors and responses.

9.5 Canonical Correlations

From the eigenvalues of $\mathbf{E}^{-1}\mathbf{H}$, the canonical correlations between the predictor and response spaces are:

$r_{c,s} = \sqrt{\frac{\lambda_s}{1 + \lambda_s}}$

The canonical $R^2_s = r_{c,s}^2$ represents the proportion of variance in the $s$-th canonical response variate explained by the corresponding canonical predictor variate. These describe the multivariate association structure.

9.6 Multivariate R² (Generalisations)

Several generalisations of $R^2$ to the multivariate case exist:

Pillai's Trace-based R²:

$R^2_{mult} = \frac{V}{s^*} = \frac{\sum_s \lambda_s/(1+\lambda_s)}{s^*}$

Wilks' Lambda-based R²:

$R^2_{Wilks} = 1 - \Lambda^{*1/t}$

Trace correlation coefficient (average $R^2$ across canonical variates):

$R^2_{trace} = \frac{\sum_s r_{c,s}^2}{s^*}$

10. Model Fit and Evaluation

10.1 Univariate R² for Each Response

Report the standard $R^2$ for each of the $q$ response regressions:

$R^2_j = 1 - \frac{\hat{\mathbf{e}}_j^T\hat{\mathbf{e}}_j}{(\mathbf{y}_j - \bar{y}_j\mathbf{1})^T(\mathbf{y}_j - \bar{y}_j\mathbf{1})}$

Adjusted R² penalises for model complexity:

$\bar{R}^2_j = 1 - (1 - R^2_j)\frac{n-1}{n-p-1}$

10.2 Trace of the Residual SSCP

The trace of the residual SSCP matrix measures total residual variation across all responses:

$\text{tr}(\hat{\mathbf{E}}^T\hat{\mathbf{E}}) = \sum_{j=1}^q \hat{\mathbf{e}}_j^T\hat{\mathbf{e}}_j = \sum_{j=1}^q SS_{res,j}$

This is the multivariate analogue of the residual sum of squares.

10.3 The Determinant Criterion ($|\hat{\boldsymbol{\Sigma}}|$)

The generalised variance $|\hat{\boldsymbol{\Sigma}}|$ (determinant of the estimated error covariance matrix) is a scalar measure of total residual variation that accounts for correlations among responses. Smaller values indicate better overall fit.

Log-determinant criterion:

$\ln|\hat{\boldsymbol{\Sigma}}| = \ln\left|\frac{\hat{\mathbf{E}}^T\hat{\mathbf{E}}}{n-p-1}\right|$

This appears in the log-likelihood and information criteria.

10.4 AIC and BIC for Multivariate Models

Akaike Information Criterion:

$AIC = n\ln|\tilde{\boldsymbol{\Sigma}}_{MLE}| + 2q(p+1)$

Or using the maximised log-likelihood:

$AIC = -2\ell(\hat{\mathbf{B}}, \tilde{\boldsymbol{\Sigma}}_{MLE}) + 2k$

Where $k = q(p+1) + q(q+1)/2$ is the total number of free parameters (coefficients plus unique elements of $\boldsymbol{\Sigma}$).

Bayesian Information Criterion:

$BIC = -2\ell(\hat{\mathbf{B}}, \tilde{\boldsymbol{\Sigma}}_{MLE}) + k\ln(n)$

Lower AIC/BIC indicates a better model relative to complexity. Use AIC for predictive accuracy; BIC for parsimony.

10.5 Model Comparison via Likelihood Ratio Test

To compare a full model (all $p$ predictors) vs. a reduced model ($p - c$ predictors, dropping the $c$ predictors specified by $\mathbf{C}$):

$\Lambda_{LR} = \frac{|\hat{\mathbf{E}}_{full}|}{|\hat{\mathbf{E}}_{reduced}|}$

$\chi^2 = -\left[n - p - 1 - \frac{m+c+1}{2}\right]\ln(\Lambda_{LR}) \sim \chi^2_{cm} \quad \text{under } H_0$

Or equivalently using Wilks' Lambda for the hypothesis being tested.

10.6 Comparing Univariate and Multivariate Fit

An important diagnostic is to compare the multivariate model fit with the fit from $q$ separate univariate regressions:

• Same coefficient estimates: $\hat{\mathbf{B}}_{MLM}$ = column-stacked OLS estimates.
• Difference in inference: MLM uses the joint error covariance structure for more powerful tests of multivariate hypotheses.
• Gain from MLM: Greatest when responses are strongly correlated and hypotheses about joint effects are of primary interest.

11. Model Diagnostics and Residuals

11.1 The Residual Matrix

The residual matrix is:

$\hat{\mathbf{E}} = \mathbf{Y} - \hat{\mathbf{Y}} = \mathbf{Y} - \mathbf{X}\hat{\mathbf{B}} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$

Each row $\hat{\boldsymbol{\epsilon}}_i = \mathbf{y}_i - \hat{\mathbf{y}}_i$ is the residual vector for observation $i$. A well-fitting model produces residuals that resemble multivariate white noise.

11.2 Univariate Residual Diagnostics (Per Response)

For each response $Y_j$, extract the $j$-th column of $\hat{\mathbf{E}}$ and apply standard univariate regression diagnostics:

Residuals vs. Fitted Values: Plot $\hat{e}_{ij}$ against $\hat{y}_{ij}$. Expect a random scatter around zero with no pattern.

Normal Q-Q Plot: Plot ordered residuals against normal quantiles. Points should fall near the diagonal.

Scale-Location Plot: Plot $\sqrt{|\hat{e}_{ij}|}$ against $\hat{y}_{ij}$. A horizontal band indicates homoscedasticity.

Residuals vs. Predictor Plots: Plot $\hat{e}_{ij}$ against each $X_k$. Patterns indicate non-linearity or omitted variable bias.

11.3 Multivariate Residual Diagnostics

11.3.1 Mahalanobis Distance of Residuals

The squared Mahalanobis distance of the residual vector for observation $i$:

$D^2_i = \hat{\boldsymbol{\epsilon}}_i^T \hat{\boldsymbol{\Sigma}}^{-1} \hat{\boldsymbol{\epsilon}}_i$

Under the model, $D^2_i / (n-p-1)$ approximately follows an $F_{q, n-p-q}$ distribution. Large $D^2_i$ indicates observation $i$ is poorly fitted on the combined response profile.

Chi-squared Q-Q plot: Plot ordered $D^2_{(i)}$ against $\chi^2_q$ quantiles. Departures from linearity indicate non-normality of multivariate residuals or outliers.

11.3.2 Standardised Multivariate Residuals

Standardise each residual vector by the estimated error covariance:

$\tilde{\boldsymbol{\epsilon}}_i = \hat{\boldsymbol{\Sigma}}^{-1/2}\hat{\boldsymbol{\epsilon}}_i$

Under correct specification, $\tilde{\boldsymbol{\epsilon}}_i \approx \mathcal{N}_q(\mathbf{0}, \mathbf{I}_q)$.

11.3.3 Cross-Response Residual Correlation

Compute the correlation matrix of the columns of $\hat{\mathbf{E}}$:

$\hat{\mathbf{R}}_\epsilon = D^{-1/2} \hat{\boldsymbol{\Sigma}} D^{-1/2}$

Where $D = \text{diag}(\hat{\boldsymbol{\Sigma}})$. This estimates the within-observation correlation structure of the errors. High correlations indicate that the responses share substantial common error variance.

11.4 Leverage and Influence

Leverage $h_{ii}$ (diagonal of $\mathbf{H}$): The same hat matrix applies to all responses in MLM. High leverage ($h_{ii} > 2(p+1)/n$) indicates an observation with unusual predictor values.

Multivariate Cook's Distance:

$D_i^{Cook} = \frac{1}{q(p+1)}\text{tr}\left[(\hat{\mathbf{B}} - \hat{\mathbf{B}}^{(-i)})^T(\mathbf{X}^T\mathbf{X})(\hat{\mathbf{B}} - \hat{\mathbf{B}}^{(-i)})\hat{\boldsymbol{\Sigma}}^{-1}\right]$

Where $\hat{\mathbf{B}}^{(-i)}$ is the coefficient matrix estimated without observation $i$. An efficient formula:

$D_i^{Cook} = \frac{h_{ii}}{q(1-h_{ii})^2} D^2_i$

Observations with $D_i^{Cook} > 1$ (or $> 4/n$) warrant investigation.

DFFITS (Multivariate):

$DFFITS_i = \frac{h_{ii}}{1-h_{ii}} D^2_i$

11.5 Testing for Multivariate Outliers

Bonferroni-corrected Mahalanobis distance test: Compare $D^2_i$ to the critical value:

$D^2_{crit} = \chi^2_{q,\, \alpha/(2n)}$

(Bonferroni-corrected at the $\alpha/(2n)$ level, e.g., $\alpha = 0.05$, $n = 100$, $q = 3$: use $\chi^2_{3, 0.00025} \approx 17.1$.)

Robust Mahalanobis distances: Use the Minimum Covariance Determinant (MCD) or Minimum Volume Ellipsoid (MVE) estimates of location and scatter, which are resistant to masking effects (outliers hiding other outliers).

11.6 Checking Multicollinearity Among Predictors

Variance Inflation Factor (VIF) for each predictor (same for all responses, since $\mathbf{X}$ is shared):

$VIF_k = \frac{1}{1 - R^2_k}$

Where $R^2_k$ is the $R^2$ from regressing $X_k$ on all other predictors.

Condition number of $\mathbf{X}^T\mathbf{X}$: Ratio of largest to smallest eigenvalue. Values $> 30$ indicate potentially problematic multicollinearity.

11.7 Checking Homoscedasticity

For each response $Y_j$, check homoscedasticity using:

• Breusch-Pagan test: Tests whether $\hat{e}_{ij}^2$ is predicted by the fitted values.
• White's test: A more general test for heteroscedasticity.
• Scale-Location plot: Visual assessment.

For the multivariate setting, Box's M test can be adapted to test whether the error covariance matrix $\boldsymbol{\Sigma}$ is the same across subgroups defined by categorical predictors.

12. Interpretation of Coefficients

12.1 Interpreting the Coefficient Matrix $\hat{\mathbf{B}}$

The coefficient matrix $\hat{\mathbf{B}}$ has dimensions $(p+1) \times q$:

• Column $j$ of $\hat{\mathbf{B}}$: The univariate regression coefficients for response $Y_j$ on all predictors. Interpretation is identical to univariate regression: $\hat{\beta}_{kj}$ is the expected change in $Y_j$ for a one-unit increase in $X_k$, holding all other predictors constant.
• Row $k$ of $\hat{\mathbf{B}}$: The multivariate coefficient vector for predictor $X_k$ — the simultaneous effect of $X_k$ on all $q$ responses. This row is tested jointly in multivariate hypothesis tests.

12.2 Ceteris Paribus Interpretation

Each coefficient $\hat{\beta}_{kj}$ represents the partial effect of $X_k$ on $Y_j$:

$\hat{\beta}_{kj} = \frac{\partial \hat{Y}_j}{\partial X_k}\Bigg|_{X_{-k} \text{ fixed}}$

This is the expected change in the mean of $Y_j$ for a one-unit increase in $X_k$, with all other predictors $X_{-k} = \{X_1, \dots, X_{k-1}, X_{k+1}, \dots, X_p\}$ held constant.

12.3 Comparing Coefficients Across Responses

A key feature of MLM is the ability to compare the effects of a predictor across multiple responses. If the responses are on the same scale (or are standardised), directly compare $\hat{\beta}_{k1}$, $\hat{\beta}_{k2}$, ..., $\hat{\beta}_{kq}$ for predictor $X_k$.

To formally test whether the effect of $X_k$ on $Y_j$ equals its effect on $Y_l$:

$H_0: \beta_{kj} = \beta_{kl} \quad \Leftrightarrow \quad H_0: \mathbf{e}_k^T\mathbf{B}\mathbf{m}_{jl} = 0$

Where $\mathbf{m}_{jl} = \mathbf{e}_j - \mathbf{e}_l$ (contrast between responses $j$ and $l$).

12.4 Interpreting Interactions

For a model with predictors $X_1$, $X_2$, and their interaction $X_1 \times X_2$:

$E[Y_j \mid X_1, X_2] = \beta_{0j} + \beta_{1j}X_1 + \beta_{2j}X_2 + \beta_{12j}X_1 X_2$

The interaction coefficient $\beta_{12j}$ represents the modification of the effect of $X_1$ on $Y_j$ per unit change in $X_2$ (and vice versa):