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

$$\begin{array}{l}MP{G}_{imk}={\beta}_{0}+{\beta}_{1}Ac{c}_{i}+{\beta}_{2}Weigh{t}_{i}+{b}_{10m}+{b}_{11}{}_{m}Ac{c}_{i}+{b}_{21k}Weigh{t}_{i}+{\epsilon}_{imk},\text{\hspace{1em}}\\ \text{\hspace{1em}}\text{\hspace{1em}}m=1,2,\mathrm{...},13,\text{\hspace{1em}}k=1,2,\mathrm{...},8,\end{array}$$

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:

$$\begin{array}{l}{b}_{1m}=\left(\begin{array}{l}{b}_{10m}\\ {b}_{11m}\end{array}\right)~N\left(0,\left(\begin{array}{cc}{\sigma}_{10}^{2}& {\sigma}_{10,11}^{}\\ {\sigma}_{10,11}^{}& {\sigma}_{11}^{2}\end{array}\right)\right),\\ b{}_{2k}\text{\hspace{0.17em}}~N\left(0,{\sigma}_{2}^{2}\right),\\ {\epsilon}_{imk}~N\left(0,{\sigma}^{2}\right).\end{array}$$

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