Class: LinearMixedModel
Extract covariance parameters of linear mixedeffects model
[
returns
the covariance parameters and related statistics in psi
,mse
,stats
]
= covarianceParameters(lme
,Name,Value
)stats
with
additional options specified by one or more Name,Value
pair
arguments.
For example, you can specify the confidence level for the confidence limits of covariance parameters.
lme
— Linear mixedeffects modelLinearMixedModel
objectLinear mixedeffects model, returned as a LinearMixedModel
object.
For properties and methods of this object, see LinearMixedModel
.
Specify optional commaseparated pairs of Name,Value
arguments.
Name
is the argument
name and Value
is the corresponding
value. Name
must appear
inside single quotes (' '
).
You can specify several name and value pair
arguments in any order as Name1,Value1,...,NameN,ValueN
.
'Alpha'
— Confidence level0.05 (default)  scalar value in the range 0 to 1Confidence level, specified as the commaseparated pair consisting
of 'Alpha'
and a scalar value in the range 0 to
1. For a value α, the confidence level is 100*(1–α)%.
For example, for 99% confidence intervals, you can specify the confidence level as follows.
Example: 'Alpha',0.01
Data Types: single
 double
psi
— Estimate of covariance parameterscell arrayEstimate of covariance parameters that parameterize the prior
covariance of the random effects, returned as a cell array of length R,
such that psi{r}
contains the covariance matrix
of random effects associated with grouping variable g_{r}, r =
1, 2, ..., R. The order of grouping variables is
the same order you enter when you fit the model.
mse
— Residual variance estimatescalar valueResidual variance estimate, returned as a scalar value.
stats
— Covariance parameter estimates and related statisticscell arrayCovariance parameter estimates and related statistics, returned as a cell array of length (R + 1) containing dataset arrays with the following columns.
Group  Grouping variable name 
Name1  Name of the first predictor variable 
Name2  Name of the second predictor variable 
Type  std (standard deviation), if Name1 and Name2 are
the same

Estimate  Standard deviation of the random effect associated with predictor Name1 or Name2 ,
if Name1 and Name2 are the sameCorrelation
between the random effects associated with predictors 
Lower  Lower limit of a 95% confidence interval for the covariance parameter 
Upper  Upper limit of a 95% confidence interval for the covariance parameter 
stats{r}
is a dataset array containing statistics
on covariance parameters for the rth grouping variable, r =
1, 2, ..., R. stats{R+1}
contains
statistics on the residual standard deviation. The dataset array for
the residual error has the fields Group
, Name
, Estimate
, Lower
,
and Upper
.
Navigate to a folder containing sample data.
cd(matlabroot)
cd('help/toolbox/stats/examples')
Load the sample data.
load fertilizer
The dataset array includes data from a splitplot experiment, where soil is divided into three blocks based on the soil type: sandy, silty, and loamy. Each block is divided into five plots, where five different types of tomato plants (cherry, heirloom, grape, vine, and plum) are randomly assigned to these plots. The tomato plants in the plots are then divided into subplots, where each subplot is treated by one of four fertilizers. This is simulated data.
Store the data in a dataset array called ds
,
for practical purposes, and define Tomato
, Soil
,
and Fertilizer
as categorical variables.
ds = fertilizer; ds.Tomato = nominal(ds.Tomato); ds.Soil = nominal(ds.Soil); ds.Fertilizer = nominal(ds.Fertilizer);
Fit a linear mixedeffects model, where Fertilizer
is
the fixedeffects variable, and the mean yield varies by the block
(soil type), and the plots within blocks (tomato types within soil
types) independently. This model corresponds to
$${y}_{ijk}={\beta}_{0}+{\displaystyle \sum _{j=2}^{5}{\beta}_{2j}I{\left[T\right]}_{ij}}+{b}_{0k}{S}_{k}+{b}_{0jk}{(S*T)}_{jk}+{\epsilon}_{ijk},$$
where i = 1, 2, ..., 60 corresponds to the observations, j = 2, ..., 5 corresponds to the tomato types, and k = 1, 2, 3 corresponds to the blocks (soil). S_{k} represents the k th soil type, and (S*T)_{jk } represents the j th tomato type nested in the k th soil type. I[T]_{ij} is the dummy variable representing the level j of the tomato type.
The random effects and observation error have the following prior distributions: b_{0k}~N(0,σ^{2}_{S}), b_{0jk}~N(0,σ^{2}_{S*T}), and ε_{ijk} ~ N(0,σ^{2}).
lme = fitlme(ds,'Yield ~ Fertilizer + (1Soil) + (1Soil:Tomato)');
Compute the covariance parameter estimates (estimates of σ^{2}_{S} and σ^{2}_{S*T}) of the randomeffects terms.
psi = covarianceParameters(lme)
psi = [4.4026e17] [ 352.8481]
Compute the residual variance (σ^{2}).
[~,mse] = covarianceParameters(lme)
mse = 151.9007
Navigate to a folder containing sample data.
cd(matlabroot)
cd('help/toolbox/stats/examples')
Load the sample data.
load weight
weight
contains data from a longitudinal
study, where 20 subjects are randomly assigned to 4 exercise programs,
and their weight loss is recorded over six 2week time periods. This
is simulated data.
Store the data in a dataset array. Define Subject
and Program
as
categorical variables.
ds = dataset(InitialWeight,Program,Subject,Week,y); ds.Subject = nominal(ds.Subject); ds.Program = nominal(ds.Program);
Fit a linear mixedeffects model where the initial weight, type of program, week, and the interaction between the week and type of program are the fixed effects. The intercept and week vary by subject.
For 'reference'
dummy variable coding, fitlme
uses
Program A as reference and creates the necessary dummy variables I[.].
This model corresponds to
$$\begin{array}{l}{y}_{im}={\beta}_{0}+{\beta}_{1}I{W}_{i}+{\beta}_{2}Wee{k}_{i}+{\beta}_{3}I{\left[PB\right]}_{i}+{\beta}_{4}I{\left[PC\right]}_{i}+{\beta}_{5}I{\left[PD\right]}_{i}\\ \text{\hspace{1em}}\text{\hspace{1em}}+b{}_{0m}+\text{\hspace{0.17em}}{b}_{1m}Wee{k}_{im}+{\epsilon}_{im},\end{array}$$
where i corresponds to the observation number, i = 1, 2, ...,120, and m corresponds to the subject number, m = 1, 2, ..., 20. β_{j} are the fixedeffects coefficients, j = 0, 1, ..., 8, and b_{0m} and b_{1m} are random effects. IW stands for initial weight and I[.] is a dummy variable representing a type of program. For example, I[PB]_{i} is the dummy variable representing Program B.
The random effects and observation error have the following prior distributions:$$\left(\begin{array}{l}{b}_{0m}\\ {b}_{1}m\end{array}\right)~N\left(0,\left(\begin{array}{cc}{\sigma}_{0}^{2}& {\sigma}_{0,1}\\ {\sigma}_{0,1}& {\sigma}_{1}^{2}\end{array}\right)\right)$$ and ε_{im} ~ N(0,σ^{2}).
lme = fitlme(ds,'y ~ InitialWeight + Program + (WeekSubject)');
Compute the estimates of covariance parameters for the random effects.
[psi,mse,stats] = covarianceParameters(lme)
psi = [2x2 double] mse = 0.0105 stats = [3x7 classreg.regr.lmeutils.titleddataset] [1x5 dataset
mse
is the estimated residual variance. It
is the estimate for σ^{2}.
To see the covariance parameters estimates for the randomeffects
terms (σ^{2}_{0}, σ^{2}_{1},
and σ^{2}_{0,1}),
index into psi
.
psi{1}
ans = 0.0572 0.0490 0.0490 0.0624
The estimate of the variance of the random effects term for the intercept, σ^{2}_{0}, is 0.0572. The estimate of the variance of the random effects term for week, σ^{2}_{1}, is 0.0624. The estimate for the covariance of the random effects terms for the intercept and week, σ_{0,1}, is 0.0490.
stats
is a 2by1 cell array. The first
cell of stats
contains the confidence intervals
for the standard deviation of the random effects and the correlation
between the random effects for intercept and week. To display them,
index into stats
.
stats{1}
ans = Covariance Type: FullCholesky Group Name1 Name2 Type Estimate Lower Upper Subject '(Intercept)' '(Intercept)' 'std' 0.23927 0.14364 0.39854 Subject 'Week' '(Intercept)' 'corr' 0.81971 0.38662 0.95658 Subject 'Week' 'Week' 'std' 0.2497 0.18303 0.34067
The display shows the name of the grouping parameter (Group
),
the randomeffects variables (Name1
, Name2
),
the type of the covariance parameters (Type
), the
estimate (Estimate
) for each parameter, and the
95% confidence intervals for the parameters (Lower
, Upper
).
The estimates in this table are related to the estimates in psi
as
follows.
The standard deviation of the randomeffects term for intercept is 0.23927 = sqrt(0.0527). Likewise, the standard deviation of the random effects term for week is 0.2497 = sqrt(0.0624). Finally, the correlation between the randomeffects terms of intercept and week is 0.81971 = 0.0490/(0.23927*0.2497).
Note that this display also shows which covariance pattern you
use when fitting the model. In this case, the covariance pattern is FullCholesky
.
To change the covariance pattern for the randomeffects terms, you
must use the 'CovariancePattern'
namevalue pair
argument when fitting the model.
The second cell of stats
includes similar
statistics for the residual standard deviation. Display the contents
of the second cell.
stats{2}
ans = Group Name Estimate Lower Upper Error 'Res Std' 0.10261 0.087882 0.11981
The estimate for residual standard deviation is the square root
of mse
, 0.10261 = sqrt(0.0105).
Load the sample data.
load carbig
Fit a linear mixedeffects 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 randomeffects 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 randomeffects term b_{1m} represents the first random effect at level m of the first grouping variable. The randomeffects 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 randomeffects term (1).
Similarly, b_{2k} stands for the second random effectsterm at level k of the second grouping variable.
σ^{2}_{10} is the variance of the randomeffects 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 randomeffects terms for the intercept and the predictor acceleration. σ^{2}_{2} is the variance of the second randomeffects term, and σ^{2} is the residual variance.
First, prepare the design matrices for fitting the linear mixedeffects 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.7989e08] 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 2by1 cell array, and stats
is a 3by1 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 randomeffects terms
for the intercept and acceleration σ_{10,11} is
–0.8387.
Now, index into the second cell of psi
.
psi{2}
ans = 6.7989e08
The second cell of psi
contains the estimate
for the variance of the randomeffects 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 randomeffects
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 offdiagonal 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.2158e05 0.00073775
The second cell of stats
has the standard
deviation estimate and the 95% confidence limits for the standard
deviation of the randomeffects 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.6466e05 0.0010229
stats{3}
ans = Group Name Estimate Lower Upper Error 'Res Std' 3.0126 2.74 3.3123