Multivariate General Linear Model

This example shows how to set up a multivariate general linear model for estimation using mvregress.

Load sample data.

This data contains measurements on a sample of 205 auto imports from 1985.

Here, model the bivariate response of city and highway MPG (columns 14 and 15).

For predictors, use wheel base (column 3), curb weight (column 7), and fuel type (column 18). The first two predictors are continuous, and for this example are centered and scaled. Fuel type is a categorical variable with two categories (11 and 20), so a dummy indicator variable is needed for the regression.

load('imports-85')
Y = X(:,14:15);
[n,d] = size(Y);

X1 = zscore(X(:,3));
X2 = zscore(X(:,7));
X3 = X(:,18)==20;

Xmat = [ones(n,1) X1 X2 X3];

The variable X3 is coded to have value 1 for the fuel type 20, and value 0 otherwise.

For convenience, the three predictors (wheel base, curb weight, and fuel type indicator) are combined into one design matrix, with an added intercept term.

Set up design matrices.

Given these predictors, the multivariate general linear model for the bivariate MPG response is

$[\begin{array}{cc} y_{11} & y_{12} \\ y_{21} & y_{22} \\ ⋮ & ⋮ \\ y_{n 1} & y_{n 2} \end{array}] = [\begin{array}{cccc} 1 & x_{11} & x_{12} & x_{13} \\ 1 & x_{21} & x_{22} & x_{23} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ 1 & x_{n 1} & x_{n 2} & x_{n 3} \end{array}] [\begin{array}{cc} β_{01} & β_{02} \\ β_{11} & β_{12} \\ β_{21} & β_{22} \\ β_{31} & β_{32} \end{array}] + [\begin{array}{cc} ϵ_{11} & ϵ_{12} \\ ϵ_{21} & ϵ_{22} \\ ⋮ & ⋮ \\ ϵ_{n 1} & ϵ_{n 2} \end{array}],$

where $ϵ_{i} = {(ϵ_{i 1}, ϵ_{i 2})}^{'} - M V N (0, Σ)$ . There are $K = 8$ regression coefficients in total.

Create a length $n = 205$ cell array of 2-by-8 (d-by-K) matrices for use with mvregress. The ith matrix in the cell array is

$X (i) = [\begin{array}{cccccccc} 1 & 0 & x_{i 1} & 0 & x_{i 2} & 0 & x_{i 3} & 0 \\ 0 & 1 & 0 & x_{i 1} & 0 & x_{i 2} & 0 & x_{i 3} \end{array}] .$

Xcell = cell(1,n);
for i = 1:n
    Xcell{i} = [kron([Xmat(i,:)],eye(d))];
end

Given this specification of the design matrices, the corresponding parameter vector is

$β = [\begin{array}{c} β_{01} \\ β_{02} \\ β_{11} \\ β_{12} \\ β_{21} \\ β_{22} \\ β_{31} \\ β_{32} \end{array}] .$

Estimate regression coefficients.

Fit the model using maximum likelihood estimation.

[beta,sigma,E,V] = mvregress(Xcell,Y);
beta

These coefficient estimates show:

The expected city and highway MPG for cars of average wheel base, curb weight, and fuel type 11 are 33.5 and 38.6, respectively. For fuel type 20, the expected city and highway MPG are 33.5476 - 9.2284 = 24.3192 and 38.5720 - 8.6663 = 29.9057.
An increase of one standard deviation in curb weight has almost the same effect on expected city and highway MPG. Given all else is equal, the expected MPG decreases by about 6.3 with each one standard deviation increase in curb weight, for both city and highway MPG.
For each one standard deviation increase in wheel base, the expected city MPG increases 0.972, while the expected highway MPG increases by only 0.395, given all else is equal.

Compute standard errors.

The standard errors for the regression coefficients are the square root of the diagonal of the variance-covariance matrix, V.

se = sqrt(diag(V))

Reshape coefficient matrix.

You can easily reshape the regression coefficients into the original 4-by-2 matrix.

B = reshape(beta,2,4)'

B = 4×2

   33.5476   38.5720
    0.9723    0.3950
   -6.3064   -6.3584
   -9.2284   -8.6663

Check model assumptions.

Under the model assumptions, $z = E Σ^{- 1 / 2}$ should be independent, with a bivariate standard normal distribution. In this 2-D case, you can assess the validity of this assumption using a scatter plot.

z = E/chol(sigma);
figure()
plot(z(:,1),z(:,2),'.')
title('Standardized Residuals')
hold on

% Overlay standard normal contours
z1 = linspace(-5,5);
z2 = linspace(-5,5);
[zx,zy] = meshgrid(z1,z2);
zgrid = [reshape(zx,100^2,1),reshape(zy,100^2,1)];
zn = reshape(mvnpdf(zgrid),100,100);
[c,h] = contour(zx,zy,zn);
clabel(c,h)

Figure contains an axes object. The axes object with title Standardized Residuals contains 2 objects of type line, contour. One or more of the lines displays its values using only markers

Several residuals are larger than expected, but overall, there is little evidence against the multivariate normality assumption.