Documentation

This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English verison of the page.

Note: This page has been translated by MathWorks. Please click here
To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

covarianceParameters

Class: LinearMixedModel

Extract covariance parameters of linear mixed-effects model

Syntax

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

Description

example

psi = covarianceParameters(lme) returns the estimated covariance parameters that parameterize the prior covariance of random effects.

example

[psi,mse] = covarianceParameters(lme) also returns an estimate of the residual variance.

example

[psi,mse,stats] = covarianceParameters(lme) also returns a cell array, stats, containing the covariance parameters and related statistics.

example

[psi,mse,stats] = covarianceParameters(lme,Name,Value) returns the covariance parameters and related statistics in 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.

Input Arguments

expand all

Linear mixed-effects model, returned as a LinearMixedModel object.

For properties and methods of this object, see LinearMixedModel.

Name-Value Pair Arguments

Specify optional comma-separated 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.

expand all

Confidence level, specified as the comma-separated 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

expand all

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 gr, r = 1, 2, ..., R. The order of grouping variables is the same order you enter when you fit the model.

Residual variance estimate, returned as a scalar value.

Covariance parameter estimates and related statistics, returned as a cell array of length (R + 1) containing dataset arrays with the following columns.

GroupGrouping variable name
Name1Name of the first predictor variable
Name2Name of the second predictor variable
Type std (standard deviation), if Name1 and Name2 are the same

corr (correlation), if Name1 and Name2 are different

EstimateStandard deviation of the random effect associated with predictor Name1 or Name2, if Name1 and Name2 are the same

Correlation between the random effects associated with predictors Name1 and Name2, if Name1 and Name2 are different

LowerLower limit of a 95% confidence interval for the covariance parameter
UpperUpper 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.

Examples

expand all

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 split-plot 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 mixed-effects model, where Fertilizer is the fixed-effects 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} + \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,\sigma^{2}_{S})$, $b_{0jk}~N(0,\sigma^{2}_{S*T}$, and $\epsilon_{ijk} ~ N(0,\sigma^{2})$.

lme = fitlme(ds,'Yield ~ Fertilizer + (1|Soil) + (1|Soil:Tomato)');

Compute the covariance parameter estimates (estimates of $\sigma^{2}_{S}$ and $\sigma^{2}_{S*T}$) of the random-effects terms.

psi = covarianceParameters(lme)
psi =

  2×1 cell array

    [2.4628e-17]
    [  352.8481]

Compute the residual variance ( $\sigma^{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 2-week 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 mixed-effects 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}{rl} y_{im} &= \beta_0 + \beta_1 IW_i + \beta_2
\mbox{Week}_i + \beta_3 I [PB]_I + \beta_4 I [PC]_i\\ &+ \beta_5 I [PD]_i
+ b_{0m} + b_{1m}\mbox{Week}_{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$. $\beta_{j}$ are the fixed-effects 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
\left( 0, \left( \begin{array}{cc} \sigma_0^2 & \sigma_{0,1}\\
\sigma_{0,1} & \sigma_1^2 \end{array} \right) \right)$$

and

$$\epsilon_{im} \sim N(0,\sigma^{2}).$$

lme = fitlme(ds,'y ~ InitialWeight + Program + (Week|Subject)');

Compute the estimates of covariance parameters for the random effects.

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

  cell

    [2×2 double]


mse =

    0.0105


stats =

  2×1 cell array

    [3×7 classreg.regr.lmeutils.titleddataset]
    [1×5 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 random-effects terms ( $\sigma^{2}_{0}$, $\sigma^{2}_{1}$, and $\sigma^{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, $\sigma^{2}_{0}$, is 0.0572. The estimate of the variance of the random effects term for week, $\sigma^{2}_{1}$, 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 2-by-1 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
    Subject    '(Intercept)'        '(Intercept)'        'std'         0.23927 
    Subject    'Week'               '(Intercept)'        'corr'        0.81971 
    Subject    'Week'               'Week'               'std'          0.2497 


    Lower      Upper  
    0.14365    0.39854
    0.38662    0.95658
    0.18303    0.34067

The display shows the name of the grouping parameter (Group), the random-effects 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 random-effects 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 random-effects 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 random-effects terms, you must use the 'CovariancePattern' name-value 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 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

$$\mbox{MPG}_{imk} = \beta_0 + \beta_1 \mbox{Acc}_i + \beta_2
\mbox{Weight}_i + b_{10m} + b_{11m} \mbox{Acc}_i + b_{21k}
\mbox{Weight}_i + \epsilon_{imk}$$

where $m = 1,2,...,13$ represents the levels for the variable Model_Year, and $k=1,2,...,8$ represents the levels for the variable Origin. $MPG_{imk}$ is the miles per gallon for the ith observation,|m| th model year, and|k| th origin that correspond to the ith observation. The random-effects terms and the observation error have the following prior distributions:

$$b_{1m} = \left( \begin{array}{c} b_{10m}\\b_{11m} \end{array} \right)
\sim 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} \sim N \left( 0, \sigma_2^2 \right),$$

$$\epsilon_{imk} \sim N \left( 0,\sigma^2 \right).$$

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.

$\sigma^{2}_{10}$ is the variance of the random-effects term for the intercept, $\sigma^{2}_{11}$ is the variance of the random effects term for the predictor acceleration, and $\sigma_{10,11}$ is the covariance of the random-effects terms for the intercept and the predictor acceleration. $\sigma^{2}_{2}$ is the variance of the second random-effects term, and $\sigma^{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 =

  2×1 cell array

    [2×2 double]
    [6.7579e-08]


mse =

    9.0802


stats =

  3×1 cell array

    [3×7 classreg.regr.lmeutils.titleddataset]
    [1×7 classreg.regr.lmeutils.titleddataset]
    [1×5 classreg.regr.lmeutils.titleddataset]

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.3081   -0.8359
   -0.8359    0.1093

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

Now, index into the second cell of psi.

psi{2}
ans =

   6.7579e-08

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

Index into the first cell of stats.

stats{1}
ans = 


    COVARIANCE TYPE: FULLCHOLESKY

    Group         Name1                 Name2                 Type      
    Model_Year    'Intercept'           'Intercept'           'std'     
    Model_Year    'Acceleration'        'Intercept'           'corr'    
    Model_Year    'Acceleration'        'Acceleration'        'std'     


    Estimate    Lower       Upper   
      2.8824      1.1116      7.4739
    -0.87704    -0.98411    -0.30222
     0.33066     0.18918     0.57796

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  
    Origin    'Weight'        'Weight'        'std'        0.00025996


    Lower         Upper    
    9.1721e-05    0.0007368

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.0133      2.8034    3.239

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      
    Model_Year    'Intercept'           'Intercept'           'std'     
    Model_Year    'Acceleration'        'Intercept'           'corr'    
    Model_Year    'Acceleration'        'Acceleration'        'std'     


    Estimate    Lower       Upper   
      2.8824       0.824      10.083
    -0.87704    -0.99176    0.018241
     0.33066     0.15873     0.68881

stats{2}
ans = 


    COVARIANCE TYPE: FULLCHOLESKY

    Group     Name1           Name2           Type         Estimate  
    Origin    'Weight'        'Weight'        'std'        0.00025996


    Lower         Upper    
    6.6115e-05    0.0010221

stats{3}
ans = 

    Group    Name             Estimate    Lower     Upper 
    Error    'Res Std'        3.0133      2.7405    3.3134

Was this topic helpful?