Load the sample data.

load carbig

Fit a linear mixed-effects model for miles per gallon
(MPG), with fixed effects for acceleration and weight, a potentially
correlated random effect for intercept and acceleration grouped by
model year, and an independent random effect for weight, grouped by
the origin of the car. This model corresponds to

where *m* represents the levels
for the variable `Model_Year`, and *k* represents
the levels for the variable `Origin`. *MPG*_{imk} is
the miles per gallon for the *i*th observation, *m*th
model year, and *k*th origin that correspond to the *i*th
observation. The random-effects terms and the observation error have
the following prior distributions:

Here, the random-effects
term *b*_{1m} represents
the first random effect at level *m* of the first
grouping variable. The random-effects term *b*_{10m} corresponds
to the first random effects term (1), for the intercept (0), at the *m*th
level (*m*) of the first grouping variable. Likewise *b*_{11m} is
the level *m* for the first predictor (1) in the
first random-effects term (1).

Similarly, *b*_{2k} stands
for the second random effects-term at level *k* of
the second grouping variable.

σ^{2}_{10} is
the variance of the random-effects term for the intercept, σ^{2}_{11} is
the variance of the random effects term for the predictor acceleration,
and σ_{10,11} is the covariance of the random-effects
terms for the intercept and the predictor acceleration. σ^{2}_{2} is
the variance of the second random-effects term, and σ^{2} is
the residual variance.

First, prepare the design matrices for fitting the linear mixed-effects
model.

X = [ones(406,1) Acceleration Weight];
Z = {[ones(406,1) Acceleration],[Weight]};
Model_Year = nominal(Model_Year);
Origin = nominal(Origin);
G = {Model_Year,Origin};

Fit the model using the design matrices.

lme = fitlmematrix(X,MPG,Z,G,'FixedEffectPredictors',....
{'Intercept','Acceleration','Weight'},'RandomEffectPredictors',...
{{'Intercept','Acceleration'},{'Weight'}},'RandomEffectGroups',{'Model_Year','Origin'});

Compute the estimates of covariance parameters for the
random effects.

[psi,mse,stats] = covarianceParameters(lme)

psi =
[2x2 double]
[6.7989e-08]
mse =
9.0755
stats =
[3x7 classreg.regr.lmeutils.titleddataset]
[1x7 classreg.regr.lmeutils.titleddataset]
[1x5 dataset ]

The residual variance `mse` is 9.0755. `psi` is
a 2-by-1 cell array, and `stats` is a 3-by-1 cell
array. To see the contents, you must index into these cell arrays.

First, index into the first cell of `psi`.

psi{1}

ans =
8.5160 -0.8387
-0.8387 0.1087

The first cell of `psi` contains the covariance
parameters for the correlated random effects for intercept σ^{2}_{10} as
8.5160, and for acceleration σ^{2}_{11} as
0.1087. The estimate for the covariance of the random-effects terms
for the intercept and acceleration σ_{10,11} is
–0.8387.

Now, index into the second cell of `psi`.

psi{2}

ans =
6.7989e-08

The second cell of `psi` contains the estimate
for the variance of the random-effects term for weight σ^{2}_{2}.

Index into the first cell of `stats`.

stats{1}

ans =
Covariance Type: FullCholesky
Group Name1 Name2 Type Estimate Lower Upper
Model_Year 'Intercept' 'Intercept' 'std' 2.9182 1.1552 7.3716
Model_Year 'Acceleration' 'Intercept' 'corr' -0.87172 -0.98267 -0.30082
Model_Year 'Acceleration' 'Acceleration' 'std' 0.32968 0.18863 0.57619

This table shows the standard deviation estimates for the random-effects
terms for intercept and acceleration. Note that the standard deviations
estimates are the square roots of the diagonal elements in the first
cell of `psi`. Specifically, 2.9182 = sqrt(8.5160)
and 0.32968 = sqrt(0.1087). The correlation is a function of the covariance
of intercept and acceleration, and the standard deviations of intercept
and acceleration. The covariance of intercept and acceleration is
the off-diagonal value in the first cell of psi, –0.8387. So,
the correlation is –.8387/(0.32968*2.92182) = –0.87.

The grouping variable for intercept and acceleration is `Model_Year`.

Index into the second cell of `stats`.

stats{2}

ans =
Covariance Type: FullCholesky
Group Name1 Name2 Type Estimate Lower Upper
Origin 'Weight' 'Weight' 'std' 0.00026075 9.2158e-05 0.00073775

The second cell of `stats` has the standard
deviation estimate and the 95% confidence limits for the standard
deviation of the random-effects term for `Weight`.
The grouping variable is `Origin`.

Index into the third cell of `stats`.

stats{3}

ans =
Group Name Estimate Lower Upper
Error 'Res Std' 3.0126 2.8028 3.238

The third cell of `stats` contains the estimate
for residual standard deviation and the 95% confidence limits. The
estimate for residual standard deviation is the square root of `mse`,
sqrt(9.0755) = 3.0126.

Construct 99% confidence intervals for the covariance
parameters.

[~,~,stats] = covarianceParameters(lme,'Alpha',0.01);
stats{1}

ans =
Covariance Type: FullCholesky
Group Name1 Name2 Type Estimate Lower Upper
Model_Year 'Intercept' 'Intercept' 'std' 2.9182 0.86341 9.8633
Model_Year 'Acceleration' 'Intercept' 'corr' -0.87172 -0.99089 0.013164
Model_Year 'Acceleration' 'Acceleration' 'std' 0.32968 0.15828 0.68669

stats{2}

ans =
Covariance Type: FullCholesky
Group Name1 Name2 Type Estimate Lower Upper
Origin 'Weight' 'Weight' 'std' 0.00026075 6.6466e-05 0.0010229

stats{3}

ans =
Group Name Estimate Lower Upper
Error 'Res Std' 3.0126 2.74 3.3123