covarianceParameters

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,Alpha=alpha)

Description

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

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

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

[psi,mse,stats] = covarianceParameters(lme,Alpha=alpha) specifies the significance level for the confidence limits of covariance parameters.

Examples

Two Random-Effects Terms for Intercept

Load and display the sample data.

load fertilizer.mat;
tbl

tbl=60×4 table
      Soil          Tomato       Fertilizer    Yield
    _________    ____________    __________    _____

    {'Sandy'}    {'Plum'    }        1          104 
    {'Sandy'}    {'Plum'    }        2          136 
    {'Sandy'}    {'Plum'    }        3          158 
    {'Sandy'}    {'Plum'    }        4          174 
    {'Sandy'}    {'Cherry'  }        1           57 
    {'Sandy'}    {'Cherry'  }        2           86 
    {'Sandy'}    {'Cherry'  }        3           89 
    {'Sandy'}    {'Cherry'  }        4           98 
    {'Sandy'}    {'Heirloom'}        1           65 
    {'Sandy'}    {'Heirloom'}        2           62 
    {'Sandy'}    {'Heirloom'}        3          113 
    {'Sandy'}    {'Heirloom'}        4           84 
    {'Sandy'}    {'Grape'   }        1           54 
    {'Sandy'}    {'Grape'   }        2           86 
    {'Sandy'}    {'Grape'   }        3           89 
    {'Sandy'}    {'Grape'   }        4          115 
      ⋮

The table tbl 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.

Convert Tomato, Soil, and Fertilizer to categorical variables.

tbl.Tomato = categorical(tbl.Tomato);
tbl.Soil = categorical(tbl.Soil);
tbl.Fertilizer = categorical(tbl.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}+{ϵ}_{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 $${ϵ}_{ijk}\sim N(0,{\sigma }^2)$$.

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

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

psi = covarianceParameters(lme)

psi=2×1 cell array
    {[3.1731e-17]}
    {[  352.8481]}

Compute the residual variance ($${\sigma }^2$$).

[~,mse] = covarianceParameters(lme)

mse = 
151.9007

Potentially Correlated Random-Effects 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 2-week time periods. This is simulated data.

Store the data in a table. Define Subject and Program as categorical variables.

tbl = table(InitialWeight,Program,Subject,Week,y);
tbl.Subject = categorical(tbl.Subject);
tbl.Program = categorical(tbl.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}{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}+{ϵ}_{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(0,\left(\begin{array}{cc}{\sigma }_0^2& {\sigma }_{0,1}\\ {\sigma }_{0,1}& {\sigma }_1^2\end{array}\right))$$

and

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

lme = fitlme(tbl,"y ~ InitialWeight + Program + (Week|Subject)");

Compute the estimates of covariance parameters for the random effects.

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

psi = 1×1 cell array
    {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 }_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 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    Lower      Upper  
    Subject    {'(Intercept)'}        {'(Intercept)'}        {'std' }        0.23927     0.14365    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 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).

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 mixed-effects 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+{ϵ}_{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, respectively. The random-effects 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),$$

$${ϵ}_{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 random-effects coefficient of Acceleration at level $$m$$.

Similarly, the coefficient $$b_{2k}$$ represents the random-effects 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 random-effects coefficient for $$b_{2k}$$, and $${\sigma }^2$$ is the residual variance.

Create design matrices for fitting the linear mixed-effects 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.0741

The residual variance mse is 9.0753. psi and stats are cell arrays that contain the covariance parameter estimates and their related statistics.

Display the first cell of psi.

psi{1}

ans = 2×2

    8.2656   -0.8697
   -0.8697    0.1158

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.

Display the second cell of psi.

psi{2}

ans = 
6.6943e-08

The second cell of psi contains the variance estimate for the random-effects coefficient corresponding to Weight.

Display the first cell of stats.

stats{1}

ans = 
    Covariance Type: FullCholesky

    Group         Name1                   Name2                   Type            Estimate    Lower       Upper    
    Model_Year    {'Intercept'   }        {'Intercept'   }        {'std' }           2.875      0.7632        10.83
    Model_Year    {'Acceleration'}        {'Intercept'   }        {'corr'}        -0.88896    -0.99313    0.0013948
    Model_Year    {'Acceleration'}        {'Acceleration'}        {'std' }         0.34031     0.16208      0.71451

The output shows the standard deviation estimates and 99% confidence bounds for the random-effects 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.

Display the second cell of stats.

stats{2}

ans = 
    Covariance Type: FullCholesky

    Group     Name1             Name2             Type           Estimate      Lower        Upper    
    Origin    {'Weight'}        {'Weight'}        {'std'}        0.00025873    6.557e-05    0.0010209

The second cell of stats contains the standard deviation estimate and 99% confidence bounds for the random-effects coefficient corresponding to Weight.

Display the third cell of stats.

stats{3}

ans = 
    Group    Name               Estimate    Lower     Upper 
    Error    {'Res Std'}        3.0123      2.7394    3.3125

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 mixed-effects model

LinearMixedModel object

Linear mixed-effects model, specified as a LinearMixedModel object constructed using fitlme or fitlmematrix.

alpha — Significance level

0.05 (default) | scalar value in the range 0 to 1

Significance level, specified as 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	std (standard deviation), if Name1 and Name2 are the same corr (correlation), if Name1 and Name2 are different
Estimate	Standard 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
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