covarianceParameters
Extract covariance parameters of linear mixedeffects model
Syntax
Description
[
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.
Examples
Two RandomEffects Terms for Intercept
Load the sample data.
load('fertilizer.mat');
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}+\sum _{j=2}^{5}{\beta}_{2j}I[T{]}_{ij}+{b}_{0jk}(S*T{)}_{jk}+{\u03f5}_{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}\sim N(0,{\sigma}_{S}^{2})$$, $${b}_{0jk}\sim N(0,{\sigma}_{S*T}^{2})$$, and $${\u03f5}_{ijk}\sim N(0,{\sigma}^{2})$$.
lme = fitlme(ds,'Yield ~ Fertilizer + (1Soil) + (1Soil:Tomato)');
Compute the covariance parameter estimates (estimates of $${\sigma}_{S}^{2}$$ and $${\sigma}_{S*T}^{2}$$) of the randomeffects terms.
psi = covarianceParameters(lme)
psi=2×1 cell array
{[2.7730e17]}
{[ 352.8481]}
Compute the residual variance ($${\sigma}^{2}$$).
[~,mse] = covarianceParameters(lme)
mse = 151.9007
Potentially Correlated RandomEffects Terms
Load the sample data.
load('weight.mat');
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}{rrl}{y}_{im}& =& {\beta}_{0}+{\beta}_{1}I{W}_{i}+{\beta}_{2}{\text{Week}}_{i}+{\beta}_{3}I[PB{]}_{I}+{\beta}_{4}I[PC{]}_{i}\\ & & +{\beta}_{5}I[PD{]}_{i}+{b}_{0m}+{b}_{1m}{\text{Week}}_{im}+{\u03f5}_{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$$. $${\beta}_{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}{c}{b}_{0m}\\ {b}_{1m}\end{array}\right)\sim N(0,\left(\begin{array}{cc}{\sigma}_{0}^{2}& {\sigma}_{0,1}\\ {\sigma}_{0,1}& {\sigma}_{1}^{2}\end{array}\right))$$
and
$${\u03f5}_{im}\sim N(0,{\sigma}^{2}).$$
lme = fitlme(ds,'y ~ InitialWeight + Program + (WeekSubject)');
Compute the estimates of covariance parameters for the random effects.
[psi,mse,stats] = covarianceParameters(lme)
psi = 1x1 cell array
{2x2 double}
mse = 0.0105
stats=2×1 cell array
{3x7 classreg.regr.lmeutils.titleddataset}
{1x5 classreg.regr.lmeutils.titleddataset}
mse
is the estimated residual variance. It is the estimate for $${\sigma}^{2}$$.
To see the covariance parameters estimates for the randomeffects terms ($${\sigma}_{0}^{2}$$, $${\sigma}_{1}^{2}$$, and $${\sigma}_{0,1}$$), index into psi
.
psi{1}
ans = 2×2
0.0572 0.0490
0.0490 0.0624
The estimate of the variance of the random effects term for the intercept, $${\sigma}_{0}^{2}$$, is 0.0572. The estimate of the variance of the random effects term for week, $${\sigma}_{1}^{2}$$, is 0.0624. The estimate for the covariance of the random effects terms for the intercept and week, $${\sigma}_{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).
Two Grouping Variables
Load the sample data.
load carbig
The variables MPG
, Acceleration
, Weight
, Model_Year
, and Origin
contain data for car mileage, acceleration, weight, year of manufacture, and the country in which the car was manufactured.
Fit a linear mixedeffects model using MPG
as the response variable, and Acceleration
and Weight
as fixed effects. Include random effects for the intercept and Acceleration
, grouped by Model_Year
, and an independent random effect for Weight
, grouped by Origin
. The formula for this model is
$${z}_{imk}={\beta}_{0}+{\beta}_{1}{x}_{i}+{\beta}_{2}{\text{y}}_{i}+{b}_{10m}+{b}_{11m}{x}_{i}+{b}_{2k}{y}_{i}+{\u03f5}_{imk}$$,
where $$x$$, $\mathit{y}$, and $\mathit{z}$ represent Acceleration
, Weight
, and MPG
, respectively. The subscript $\mathit{i}$ corresponds to the row in MPG
for the observation, and the subscripts $$m=1,2,...,13$$ and $$k=1,2,...,8$$ correspond to the levels for Model_Year
and Origin
. The randomeffects coefficients and the observation error have the following prior distributions:
$${b}_{1m}=\left(\begin{array}{c}{b}_{10m}\\ {b}_{11m}\end{array}\right)\sim N(0,\left(\begin{array}{cc}{\sigma}_{10}^{2}& {\sigma}_{10,11}\\ {\sigma}_{10,11}& {\sigma}_{11}^{2}\end{array}\right)),$$
$${b}_{2k}\sim N(0,{\sigma}_{2}^{2}),$$
$${\u03f5}_{imk}\sim N(0,{\sigma}^{2}).$$
The coefficient vector $${b}_{1m}$$ represents the random effect of Model_Year
at level $$m$$.
$${b}_{10m}$$ is the random intercept at level $$m$$.
$${b}_{11m}$$ is the randomeffects coefficient of
Acceleration
at level $$m$$.
Similarly, the coefficient $${b}_{2k}$$ represents the randomeffects coefficient for Weight
at level $$k$$ of Origin
.
$${\sigma}_{10}^{2}$$ is the variance of $${b}_{10m}$$, $${\sigma}_{11}^{2}$$ is the variance of the random effects coefficient $${b}_{11m}$$, and $${\sigma}_{10,11}$$ is the covariance of $${b}_{10m}$$ and $${b}_{11m}$$. $${\sigma}_{2}^{2}$$ is the variance of the randomeffects coefficient for $${b}_{2k}$$, and $${\sigma}^{2}$$ is the residual variance.
Create design matrices for fitting the linear mixedeffects model.
F = [ones(406,1) Acceleration Weight]; R = {[ones(406,1) Acceleration],Weight}; Model_Year = nominal(Model_Year); Origin = nominal(Origin); G = {Model_Year,Origin};
Fit the model using F
as the fixed effects, MPG
as the response, R
as the random effects, and G
as the grouping variables.
lme = fitlmematrix(F,MPG,R,G,'FixedEffectPredictors',.... {'Intercept','Acceleration','Weight'},'RandomEffectPredictors',... {{'Intercept','Acceleration'},{'Weight'}},'RandomEffectGroups',{'Model_Year','Origin'});
Calculate covariance parameter estimates and 99% confidence intervals for the random effects. Display the mean squared error for the residual variance.
[psi,mse,stats] = covarianceParameters(lme,Alpha=0.01); mse
mse = 9.0753
The residual variance mse
is 9.0753
. psi
and stats
are cell arrays that contain covariance parameter estimates and their related statistics.
Inspect the first cell of psi
.
psi{1}
ans = 2×2
8.2649 0.8698
0.8698 0.1157
The first cell of psi
contains the covariance parameter estimates for the correlated random effects coefficients. The variance estimate corresponding to the intercept is 8.2649
, and the variance estimate corresponding to Acceleration
is 0.1157
. The covariance estimate corresponding to the intercept and Acceleration
is 0.8698
.
Inspect the second cell of psi
.
psi{2}
ans = 6.6770e08
The second cell of psi
contains the variance estimate for the randomeffects coefficient corresponding to Weight
.
Inspect the first cell of stats
.
stats{1}
ans = COVARIANCE TYPE: FULLCHOLESKY Group Name1 Name2 Type Estimate Lower Upper Model_Year {'Intercept' } {'Intercept' } {'std' } 2.8749 0.76378 10.821 Model_Year {'Acceleration'} {'Intercept' } {'corr'} 0.88949 0.99322 0.0026856 Model_Year {'Acceleration'} {'Acceleration'} {'std' } 0.34015 0.16213 0.71364
The output shows the standard deviation estimates and 99% confidence bounds for the randomeffects coefficients corresponding to the intercept and Acceleration
. The output also displays the name of the grouping variable, Model_Year
. Note that the standard deviation estimates are the square roots of the diagonal elements in the first cell of psi
.
Inspect the second cell of stats
.
stats{2}
ans = COVARIANCE TYPE: FULLCHOLESKY Group Name1 Name2 Type Estimate Lower Upper Origin {'Weight'} {'Weight'} {'std'} 0.0002584 6.5446e05 0.0010202
The second cell of stats
contains the standard deviation estimate and 99% confidence bounds for the random effects coefficient corresponding to Weight
.
Inspect the third cell of stats
.
stats{3}
ans = Group Name Estimate Lower Upper Error {'Res Std'} 3.0125 2.7395 3.3127
The third cell of stats
contains the residual standard deviation estimate and corresponding 99% confidence bounds. Note that the residual standard deviation estimate is the square root of mse
.
Input Arguments
lme
— Linear mixedeffects model
LinearMixedModel
object
Linear mixedeffects model, specified as a LinearMixedModel
object constructed using fitlme
or fitlmematrix
.
NameValue Arguments
Specify optional pairs of arguments as
Name1=Value1,...,NameN=ValueN
, where Name
is
the argument name and Value
is the corresponding value.
Namevalue arguments must appear after other arguments, but the order of the
pairs does not matter.
Before R2021a, use commas to separate each name and value, and enclose
Name
in quotes.
Example: [psi,mse,stats] =
covarianceParameters(lme,'Alpha',0.01);
Alpha
— Significance level
0.05 (default)  scalar value in the range 0 to 1
Significance 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
Output Arguments
psi
— Estimate of covariance parameters
cell array
Estimate 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 estimate
scalar value
Residual variance estimate, returned as a scalar value.
stats
— Covariance parameter estimates and related statistics
cell array
Covariance 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 

Estimate 
Standard deviation of the random effect associated
with predictor Correlation 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
.
Version History
Introduced in R2013b
See Also
MATLAB Command
You clicked a link that corresponds to this MATLAB command:
Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
 América Latina (Español)
 Canada (English)
 United States (English)
Europe
 Belgium (English)
 Denmark (English)
 Deutschland (Deutsch)
 España (Español)
 Finland (English)
 France (Français)
 Ireland (English)
 Italia (Italiano)
 Luxembourg (English)
 Netherlands (English)
 Norway (English)
 Österreich (Deutsch)
 Portugal (English)
 Sweden (English)
 Switzerland
 United Kingdom (English)