Accelerating the pace of engineering and science

# fitlmematrix

Fit linear mixed-effects model

## Syntax

• lme = fitlmematrix(X,y,Z,[]) example
• lme = fitlmematrix(X,y,Z,G) example
• lme = fitlmematrix(___,Name,Value) example

## Description

example

lme = fitlmematrix(X,y,Z,[]) creates a linear mixed-effects model of the responses y using the fixed-effects design matrix X and random-effects design matrix or matrices in Z.

[] implies that there is one group. That is, the grouping variable G is ones(n,1), where n is the number of observations. Using fitlmematrix(X,Y,Z,[]) without a specified covariance pattern most likely results in a nonidentifiable model. This syntax is recommended only if you build the grouping information into the random effects design Z and specify a covariance pattern for the random effects using the 'CovariancePattern' name-value pair argument.

example

lme = fitlmematrix(X,y,Z,G) creates a linear mixed-effects model of the responses y using the fixed-effects design matrix X and random-effects design matrix Z or matrices in Z, and the grouping variable or variables in G.

example

lme = fitlmematrix(___,Name,Value) also creates a linear mixed-effects model with additional options specified by one or more Name,Value pair arguments, using any of the previous input arguments.

For example, you can specify the names of the response, predictor, and grouping variables. You can also specify the covariance pattern, fitting method, or the optimization algorithm.

## Examples

expand all

### No Grouping Variable Specified

`load carsmall`

Fit a linear mixed-effects model, where miles per gallon (MPG) is the response, weight is the predictor variable, and the intercept varies by model year. First, define the design matrices. Then, fit the model using the specified design matrices.

```y = MPG;
X = [ones(size(Weight)), Weight];
Z = ones(size(y));
lme = fitlmematrix(X,y,Z,Model_Year)```
```lme =

Linear mixed-effects model fit by ML

Model information:
Number of observations              94
Fixed effects coefficients           2
Random effects coefficients          3
Covariance parameters                2

Formula:
y ~ x1 + x2 + (z11 | g1)

Model fit statistics:
AIC       BIC       LogLikelihood    Deviance
486.09    496.26    -239.04          478.09

Fixed effects coefficients (95% CIs):
Name        Estimate      SE           tStat      DF    pValue        Lower         Upper
'x1'            43.575       2.3038     18.915    92    1.8371e-33            39        48.151
'x2'        -0.0067097    0.0004242    -15.817    92    5.5373e-28    -0.0075522    -0.0058672

Random effects covariance parameters (95% CIs):
Group: g1 (3 Levels)
Name1        Name2        Type         Estimate    Lower     Upper
'z11'        'z11'        'std'        3.301       1.4448    7.5421

Group: Error
Name             Estimate    Lower     Upper
'Res Std'        2.8997      2.5075    3.3532```

Now, fit the same model by building the grouping into the Z matrix.

```Z = double([Model_Year==70, Model_Year==76, Model_Year==82]);
lme = fitlmematrix(X,y,Z,[],'Covariancepattern','Isotropic')```
```lme =

Linear mixed-effects model fit by ML

Model information:
Number of observations              94
Fixed effects coefficients           2
Random effects coefficients          3
Covariance parameters                2

Formula:
y ~ x1 + x2 + (z11 + z12 + z13 | g1)

Model fit statistics:
AIC       BIC       LogLikelihood    Deviance
486.09    496.26    -239.04          478.09

Fixed effects coefficients (95% CIs):
Name        Estimate      SE           tStat      DF    pValue        Lower         Upper
'x1'            43.575       2.3038     18.915    92    1.8371e-33            39        48.151
'x2'        -0.0067097    0.0004242    -15.817    92    5.5373e-28    -0.0075522    -0.0058672

Random effects covariance parameters (95% CIs):
Group: g1 (1 Levels)
Name1        Name2        Type         Estimate    Lower     Upper
'z11'        'z11'        'std'        3.301       1.4448    7.5421

Group: Error
Name             Estimate    Lower     Upper
'Res Std'        2.8997      2.5075    3.3532
```

### Longitudinal Study with a Covariate

Navigate to a folder containing sample data.

```cd(matlabroot)
cd('help/toolbox/stats/examples')```

`load weight`

weight contains data from a longitudinal study, where 20 subjects are randomly assigned 4 exercise programs (A, B, C, D) and their weight loss is recorded over six 2-week time periods. This is simulated data.

Define Subject and Program as categorical variables. Create the design matrices for a linear mixed-effects model, with the initial weight, type of program, week, and the interaction between the week and type of program as the fixed effects. The intercept and coefficient of week vary by subject.

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{ }\text{ }+{\beta }_{6}\left(Wee{k}_{i}*I{\left[PB\right]}_{i}\right)+{\beta }_{7}\left(Wee{k}_{i}*I{\left[PC\right]}_{i}\right)+{\beta }_{8}\left(Wee{k}_{i}*I{\left[PD\right]}_{i}\right)\\ \text{ }\text{ }+b{}_{0m}+\text{\hspace{0.17em}}{b}_{1m}Wee{k}_{im}+{\epsilon }_{im},\end{array}$

where i = 1, 2, ..., 120, and m = 1, 2, ..., 20. βj are the fixed-effects coefficients, j = 0, 1, ..., 8, and b0m and b1m 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 type B. The random effects and observation error have these prior distributions: b0m~N(0,σ20), b1m~N(0,σ21), and εim ~ N(0,σ2).

```Subject = nominal(Subject);
Program = nominal(Program);
D = dummyvar(Program); % Create dummy variables for Program
X = [ones(120,1), InitialWeight, D(:,2:4), Week,...
D(:,2).*Week, D(:,3).*Week, D(:,4).*Week];
Z = [ones(120,1), Week];
G = Subject;```

Since the model has an intercept, you only need the dummy variables for programs B, C, and D. This is also known as the 'reference' method of coding dummy variables.

Fit the model using fitlmematrix with the defined design matrices and grouping variables.

```lme = fitlmematrix(X,y,Z,G,'FixedEffectPredictors',...
{'Intercept','InitWeight','PrgB','PrgC','PrgD','Week','Week_PrgB','Week_PrgC','Week_PrgD'},...
'RandomEffectPredictors',{{'Intercept','Week'}},'RandomEffectGroups',{'Subject'})```
```lme =

Linear mixed-effects model fit by ML

Model information:
Number of observations             120
Fixed effects coefficients           9
Random effects coefficients         40
Covariance parameters                4

Formula:
Linear Mixed Formula with 10 predictors.

Model fit statistics:
AIC        BIC       LogLikelihood    Deviance
-22.981    13.257    24.49            -48.981

Fixed effects coefficients (95% CIs):
Name                Estimate     SE           tStat       DF     pValue       Lower         Upper
'Intercept'           0.66105      0.25892      2.5531    111     0.012034       0.14798       1.1741
'InitWeight'        0.0031879    0.0013814      2.3078    111     0.022863    0.00045067    0.0059252
'PrgB'                0.36079      0.13139       2.746    111    0.0070394       0.10044      0.62113
'PrgC'              -0.033263      0.13117    -0.25358    111      0.80029      -0.29319      0.22666
'PrgD'                0.11317      0.13132     0.86175    111      0.39068      -0.14706       0.3734
'Week'                 0.1732     0.067454      2.5677    111     0.011567      0.039536      0.30686
'Week_PrgB'          0.038771     0.095394     0.40644    111      0.68521      -0.15026       0.2278
'Week_PrgC'          0.030543     0.095394     0.32018    111      0.74944      -0.15849      0.21957
'Week_PrgD'          0.033114     0.095394     0.34713    111      0.72915      -0.15592      0.22214

Random effects covariance parameters (95% CIs):
Group: Subject (20 Levels)
Name1              Name2              Type          Estimate    Lower      Upper
'Intercept'        'Intercept'        'std'         0.18407     0.12281    0.27587
'Week'             'Intercept'        'corr'        0.66841     0.21076    0.88573
'Week'             'Week'             'std'         0.15033     0.11004    0.20537

Group: Error
Name             Estimate    Lower       Upper
'Res Std'        0.10261     0.087882    0.11981```

The p-values 0.0228 and 0.0115 indicate significant effects of the initial weights of the subjects and the time factor in the amount of weight lost. The weight loss of subjects who are in program B is significantly different relative to the weight loss of subjects who are in program A. The lower and upper limits of the covariance parameters for the random effects do not include zero, thus they seem significant. You can also test the significance of the random-effects using the compare method.

### Random Intercept Model

`load flu`

The flu dataset array has a Date variable, and 10 variables for estimated influenza rates (in 9 different regions, estimated from Google® searches, plus a nationwide estimate from the Centers for Disease Control and Prevention, CDC).

To fit a linear-mixed effects model, where the influenza rates are the responses, combine the nine columns corresponding to the regions into a tall array that has a single response variable, FluRate, and a nominal variable, Region, the nationwide estimate WtdILI, that shows which region each estimate is from, and the grouping variable Date.

```flu2 = stack(flu,2:10,'NewDataVarName','FluRate',...
'IndVarName','Region');
flu2.Date = nominal(flu2.Date);```

Define the design matrices for a random-intercept linear mixed-effects model, where the intercept varies by Date. The corresponding model is

${y}_{im}={\beta }_{0}+{\beta }_{1}WtdIL{I}_{im}+{b}_{0m}+{\epsilon }_{im},\text{ }i=1,2,...,468,\text{ }m=1,2,...,52,$

where yim is the observation i for level m of grouping variable Date. b0m is the random effect for level m of the grouping variable Date and εim is the observation error for observation i. The random effect has the prior distribution, b0m ~ N(0,σ2FR) and the error term has the distribution, εim ~ N(0,σ2).

```y = flu2.FluRate;
X = [ones(468,1) flu2.WtdILI];
Z = [ones(468,1)];
G = flu2.Date;```

Fit the linear mixed-effects model.

```lme = fitlmematrix(X,y,Z,G,'FixedEffectPredictors',{'Intercept','NationalRate'},...
'RandomEffectPredictors',{{'Intercept'}},'RandomEffectGroups',{'Date'})```
```lme =

Linear mixed-effects model fit by ML

Model information:
Number of observations             468
Fixed effects coefficients           2
Random effects coefficients         52
Covariance parameters                2

Formula:
y ~ Intercept + NationalRate + (Intercept | Date)

Model fit statistics:
AIC       BIC       LogLikelihood    Deviance
286.24    302.83    -139.12          278.24

Fixed effects coefficients (95% CIs):
Name                  Estimate    SE          tStat     DF     pValue        Lower       Upper
'Intercept'           0.16385     0.057525    2.8484    466     0.0045885    0.050813    0.27689
'NationalRate'         0.7236     0.032219    22.459    466    3.0502e-76     0.66028    0.78691

Random effects covariance parameters (95% CIs):
Group: Date (52 Levels)
Name1              Name2              Type         Estimate    Lower      Upper
'Intercept'        'Intercept'        'std'        0.17146     0.13227    0.22226

Group: Error
Name             Estimate    Lower      Upper
'Res Std'        0.30201     0.28217    0.32324```

The confidence limits of the standard deviation of the random-effects term σ2b0m, do not include zero (0.13227, 0.22226), which indicates that the random-effects term is significant. You can also test the significance of the random-effects using compare method.

The estimated value of an observation is the sum of the fixed-effects values and value of the random effect at the grouping variable level corresponding to that observation. For example, the estimated flu rate for observation 28

$\begin{array}{l}{\stackrel{^}{y}}_{28}={\stackrel{^}{\beta }}_{0}+{\stackrel{^}{\beta }}_{1}WtdIL{I}_{28}+{\stackrel{^}{b}}_{10/30/2005}\\ \text{ }\text{ }\text{ }\text{ }\text{ }=0.1639+0.7236*\left(1.343\right)+0.3318\\ \text{ }\text{ }\text{ }\text{ }\text{ }=1.46749,\end{array}$

where $\stackrel{^}{b}$ is the best linear unbiased predictor (BLUP) of the random effects for the intercept. You can compute this value as follows.

```beta = fixedEffects(lme);
[~,~,STATS] = randomEffects(lme); % compute the random effects statistics STATS
STATS.Level = nominal(STATS.Level);
y_hat = beta(1) + beta(2)*flu2.WtdILI(28) + STATS.Estimate(STATS.Level=='10/30/2005')```
```y_hat =

1.4674```

You can simply display the fitted value using the fitted(lme) method.

```F = fitted(lme);
F(28)```
```ans =

1.4674```

### Randomized Block Design

Navigate to a folder containing sample data.

```cd(matlabroot)
cd('help/toolbox/stats/examples')```

`load shift`

The data shows the deviations from the target quality characteristic measured from the products that five operators manufacture during three shifts: morning, evening, and night. This is a randomized block design, where the operators are the blocks. The experiment is designed to study the impact of the time of shift on the performance. The performance measure is the deviations of the quality characteristics from the target value. This is simulated data.

Define the design matrices for a linear mixed-effects model with a random intercept grouped by operator, and shift as the fixed effects. Use the 'effects' contrasts. 'effects' contrasts mean that the coefficients sum to 0. You need to create two contrast coded variables in the fixed-effects design matrix, X1 and X2, where

$Shift\text{_}Evening=\left\{\begin{array}{c}0,\text{ }if\text{\hspace{0.17em}}Morning\\ 1,\text{ }if\text{\hspace{0.17em}}Evening\\ -1,\text{ }if\text{\hspace{0.17em}}Night\text{ }\end{array}\text{ }and\text{ }Shift\text{_}Morning=\left\{\begin{array}{c}1,\text{ }if\text{\hspace{0.17em}}Morning\\ 0,\text{ }if\text{\hspace{0.17em}}Evening\\ -1,\text{ }if\text{\hspace{0.17em}}Night\text{ }\end{array}.$

The model corresponds to

where i represents the observations, and m represents the operators, i = 1, 2, ..., 15, and m = 1, 2, ..., 5. The random effects and the observation error have these distributions: b0m ~ N(0,σ2b0m) and εim ~ N(0,σ2).

```S = shift.Shift;
X1 = (S=='Morning') - (S=='Night');
X2 = (S=='Evening') - (S=='Night');
X = [ones(15,1), X1, X2];
y = shift.QCDev;
Z = ones(15,1);
G = shift.Operator;```

Fit a linear mixed-effects model using the specified design matrices and restricted maximum likelihood method.

```lme = fitlmematrix(X,y,Z,G,'FitMethod','REML','FixedEffectPredictors',....
{'Intercept','S_Morning','S_Evening'},'RandomEffectPredictors',{{'Intercept'}},...
'RandomEffectGroups',{'Operator'},'DummyVarCoding','effects')```
```lme =

Linear mixed-effects model fit by REML

Model information:
Number of observations              15
Fixed effects coefficients           3
Random effects coefficients          5
Covariance parameters                2

Formula:
y ~ Intercept + S_Morning + S_Evening + (Intercept | Operator)

Model fit statistics:
AIC       BIC       LogLikelihood    Deviance
58.913    61.337    -24.456          48.913

Fixed effects coefficients (95% CIs):
Name               Estimate    SE         tStat      DF    pValue       Lower      Upper
'Intercept'          3.6525    0.94109     3.8812    12    0.0021832     1.6021       5.703
'S_Morning'        -0.91973    0.31206    -2.9473    12     0.012206    -1.5997    -0.23981
'S_Evening'        -0.53293    0.31206    -1.7078    12      0.11339    -1.2129     0.14699

Random effects covariance parameters (95% CIs):
Group: Operator (5 Levels)
Name1              Name2              Type         Estimate    Lower      Upper
'Intercept'        'Intercept'        'std'        2.0457      0.98207    4.2612

Group: Error
Name             Estimate    Lower      Upper
'Res Std'        0.85462     0.52357    1.395
```

Compute the best linear unbiased predictor (BLUP) estimates of random effects.

`B = randomEffects(lme)`
```B =

0.5775
1.1757
-2.1715
2.3655
-1.9472```

The estimated deviation from the target quality characteristics for the third operator working the evening shift is

$\begin{array}{l}{\stackrel{^}{y}}_{Evening,Operator3}={\stackrel{^}{\beta }}_{0}+{\stackrel{^}{\beta }}_{1}Shift\text{_}Evening+{\stackrel{^}{b}}_{03}\\ \text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{\hspace{0.17em}}=3.6525-0.53293-2.1715\\ \text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{\hspace{0.17em}}=0.94807.\end{array}$

You can also display this value as follows.

```F = fitted(lme);
F(shift.Shift=='Evening' & shift.Operator=='3')```
```ans =

0.9481```

### Correlated and Uncorrelated Random-Effects Terms

`load carbig`

Fit a linear mixed-effects model for miles per gallon (MPG), with fixed effects for acceleration and horsepower, and uncorrelated random effect for intercept and acceleration grouped by the model year. This model corresponds to

$MP{G}_{im}={\beta }_{0}+{\beta }_{1}Ac{c}_{i}+{\beta }_{2}HP+{b}_{0m}+{b}_{1}{}_{m}Ac{c}_{im}+{\epsilon }_{im},\text{ }m=1,2,3,$

with the random-effects terms having these distributions: b0m ~ N(0,σ20), and b1m ~ N(0,σ21). m represents the model year.

First, prepare the design matrices for fitting the linear mixed-effects model.

```X = [ones(406,1) Acceleration Horsepower];
Z = {ones(406,1),Acceleration};
G = {Model_Year,Model_Year};
Model_Year = nominal(Model_Year);
```

Now, fit the model using fitlmematrix with the defined design matrices and grouping variables.

```lme = fitlmematrix(X,MPG,Z,G,'FixedEffectPredictors',....
{'Intercept','Acceleration','Horsepower'},'RandomEffectPredictors',...
{{'Intercept'},{'Acceleration'}},'RandomEffectGroups',{'Model_Year','Model_Year'})
```
```lme =

Linear mixed-effects model fit by ML

Model information:
Number of observations             392
Fixed effects coefficients           3
Random effects coefficients         26
Covariance parameters                3

Formula:
y ~ Intercept + Acceleration + Horsepower + (Intercept | Model_Year) + (Acceleration | Model_Year)

Model fit statistics:
AIC       BIC       LogLikelihood    Deviance
2194.5    2218.3    -1091.3          2182.5

Fixed effects coefficients (95% CIs):
Name                  Estimate    SE           tStat      DF     pValue        Lower       Upper
'Intercept'             49.839       2.0518     24.291    389    5.6168e-80      45.806      53.873
'Acceleration'        -0.58565      0.10846    -5.3995    389    1.1652e-07     -0.7989     -0.3724
'Horsepower'          -0.16534    0.0071227    -23.213    389    1.9755e-75    -0.17934    -0.15133

Random effects covariance parameters (95% CIs):
Group: Model_Year (13 Levels)
Name1              Name2              Type         Estimate      Lower    Upper
'Intercept'        'Intercept'        'std'        8.0669e-07    NaN      NaN

Group: Model_Year (13 Levels)
Name1                 Name2                 Type         Estimate    Lower      Upper
'Acceleration'        'Acceleration'        'std'        0.18783     0.12523    0.28172

Group: Error
Name             Estimate    Lower     Upper
'Res Std'        3.7258      3.4698    4.0007```

Note that the random effects covariance parameters for intercept and acceleration are separate in the display. The standard deviation of the random effect for the intercept does not seem significant.

Refit the model with potentially correlated random effects for intercept and acceleration. In this case, the random-effects terms has this prior distribution

${b}_{m}=\left(\begin{array}{l}{b}_{0m}\\ {b}_{1m}\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),$

where m represents the model year.

First, prepare the random-effects design matrix and grouping variable.

```Z = [ones(406,1) Acceleration];
G = Model_Year;

lme = fitlmematrix(X,MPG,Z,G,'FixedEffectPredictors',....
{'Intercept','Acceleration','Horsepower'},'RandomEffectPredictors',...
{{'Intercept','Acceleration'}},'RandomEffectGroups',{'Model_Year'})
```
```lme =

Linear mixed-effects model fit by ML

Model information:
Number of observations             392
Fixed effects coefficients           3
Random effects coefficients         26
Covariance parameters                4

Formula:
y ~ Intercept + Acceleration + Horsepower + (Intercept + Acceleration | Model_Year)

Model fit statistics:
AIC       BIC       LogLikelihood    Deviance
2193.5    2221.3    -1089.7          2179.5

Fixed effects coefficients (95% CIs):
Name                  Estimate    SE           tStat      DF     pValue        Lower       Upper
'Intercept'             50.133       2.2652     22.132    389    7.7727e-71      45.679      54.586
'Acceleration'        -0.58327      0.13394    -4.3545    389    1.7075e-05    -0.84661    -0.31992
'Horsepower'          -0.16954    0.0072609     -23.35    389     5.188e-76    -0.18382    -0.15527

Random effects covariance parameters (95% CIs):
Group: Model_Year (13 Levels)
Name1                 Name2                 Type          Estimate    Lower       Upper
'Intercept'           'Intercept'           'std'           3.3475      1.2862      8.7119
'Acceleration'        'Intercept'           'corr'        -0.87971    -0.98501    -0.29675
'Acceleration'        'Acceleration'        'std'          0.33789      0.1825     0.62558

Group: Error
Name             Estimate    Lower     Upper
'Res Std'        3.6874      3.4298    3.9644
```

Note that the random effects covariance parameters for intercept and acceleration are together in the display, with an addition of the correlation between the intercept and acceleration. The confidence intervals for the standard deviations and the correlation between the random effects for intercept and acceleration do not include 0s, hence they seem significant. You can compare these two models using the compare method.

### Specify the Covariance Pattern

Navigate to a folder containing sample data.

```cd(matlabroot)
cd('help/toolbox/stats/examples')```

`load weight`

weight contains data from a longitudinal study, where 20 subjects are randomly assigned 4 exercise programs, and their weight loss is recorded over six 2-week time periods. This is simulated data.

Define Subject and Program as categorical variables.

```Subject = nominal(Subject);
Program = nominal(Program);
```

Create the design matrices for a linear mixed-effects model, with the initial weight, type of program, and week as the fixed effects.

```D = dummyvar(Program);
X = [ones(120,1), InitialWeight, D(:,2:4), Week];
Z = [ones(120,1) Week];
G = Subject;```

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{ }\text{ }+b{}_{0m}+\text{\hspace{0.17em}}{b}_{1m}Week{2}_{im}+{b}_{2m}Week{4}_{im}+{b}_{3m}Week{6}_{im}+{b}_{4m}Week{8}_{im}\\ \text{ }\text{ }+{b}_{5m}Week{10}_{im}+{b}_{6m}Week{12}_{im}+{\epsilon }_{im},\end{array}$

where i = 1, 2, ..., 120, and m = 1, 2, ..., 20.

βj are the fixed-effects coefficients, j = 0, 1, ..., 8, and b1m and b1m 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 type B. The random effects and observation error have these prior distributions: b0m~N(0,σ20), b1m~N(0,σ21), and εim ~ N(0,σ2).

Fit the model using fitlmematrix with the defined design matrices and grouping variables. Assume the repeated observations collected on a subject have common variance along diagonals.

```lme = fitlmematrix(X,y,Z,G,'FixedEffectPredictors',...
{'Intercept','InitWeight','PrgB','PrgC','PrgD','Week'},...
'RandomEffectPredictors',{{'Intercept','Week'}},...
'RandomEffectGroups',{'Subject'},'CovariancePattern','Isotropic')```
```lme =

Linear mixed-effects model fit by ML

Model information:
Number of observations             120
Fixed effects coefficients           6
Random effects coefficients         40
Covariance parameters                2

Formula:
y ~ Intercept + InitWeight + PrgB + PrgC + PrgD + Week + (Intercept + Week | Subject)

Model fit statistics:
AIC        BIC       LogLikelihood    Deviance
-24.783    -2.483    20.391           -40.783

Fixed effects coefficients (95% CIs):
Name                Estimate     SE           tStat       DF     pValue        Lower        Upper
'Intercept'            0.4208      0.28169      1.4938    114       0.13799     -0.13723      0.97883
'InitWeight'        0.0045552    0.0015338      2.9699    114     0.0036324    0.0015168    0.0075935
'PrgB'                0.36993      0.12119      3.0525    114     0.0028242      0.12986         0.61
'PrgC'              -0.034009       0.1209    -0.28129    114       0.77899     -0.27351       0.2055
'PrgD'                  0.121      0.12111     0.99911    114       0.31986     -0.11891      0.36091
'Week'                0.19881     0.037134      5.3538    114    4.5191e-07      0.12525      0.27237

Random effects covariance parameters (95% CIs):
Group: Subject (20 Levels)
Name1              Name2              Type         Estimate    Lower      Upper
'Intercept'        'Intercept'        'std'        0.16561     0.12896    0.21269

Group: Error
Name             Estimate    Lower       Upper
'Res Std'        0.10272     0.088014    0.11987 ```

## Input Arguments

expand all

### X — Fixed-effects design matrixn-by-p matrix

Fixed-effects design matrix, specified as an n-by-p matrix, where n is the number of observations, and p is the number of fixed-effects predictor variables. Each row of X corresponds to one observation, and each column of X corresponds to one variable.

Data Types: single | double

### y — Response valuesn-by-1 vector

Response values, specified as an n-by-1 vector, where n is the number of observations.

Data Types: single | double

### Z — Random-effects designn-by-q matrix | cell array of R n-by-q(r) matrices, r = 1, 2, ..., R

Random-effects design, specified as either of the following.

• If there is one random-effects term in the model, then Z must be an n-by-q matrix, where n is the number of observations and q is the number of variables in the random-effects term.

• If there are R random-effects terms, then Z must be a cell array of length R. Each cell of Z contains an n-by-q(r) design matrix Z{r}, r = 1, 2, ..., R, corresponding to each random-effects term. Here, q(r) is the number of random effects term in the rth random effects design matrix, Z{r}.

Data Types: single | double | cell

### G — Grouping variable or variablesn-by-1 vector | cell array of R n-by-1 vectors

Grouping variable or variables, specified as either of the following.

• If there is one random-effects term, then G must be an n-by-1 vector corresponding to a single grouping variable with M levels or groups.

G can be a categorical vector, numeric vector, character array, or cell array of strings.

• If there are multiple random-effects terms, then G must be a cell array of length R. Each cell of G contains a grouping variable G{r}, r = 1, 2, ..., R, with M(r) levels.

G{r} can be a categorical vector, numeric vector, character array, or cell array of strings.

Data Types: single | double | char | cell

### 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.

Example: 'CovariancePattern','Diagonal','DummyVarCoding','full','Optimizer','fminunc' specifies a random-effects covariance pattern with zero off-diagonal elements, creates a dummy variable for each level of a categorical variable, and uses the fminunc optimization algorithm.

### 'FixedEffectPredictors' — Names of columns in fixed-effects design matrix{'x1','x2',...,'xP'} (default) | cell array of length p

Names of columns in the fixed-effects design matrix X, specified as the comma-separated pair consisting of 'FixedEffectPredictors' and a cell array of length p.

For example, if you have a constant term and two predictors, say TimeSpent and Gender, where Female is the reference level for Gender, as the fixed effects, then you can specify the names of your fixed effects in the following way. Gender_Male represents the dummy variable you must create for category Male. You can choose different names for these variables.

Example: 'FixedEffectPredictors',{'Intercept','TimeSpent','Gender_Male'},

Data Types: cell

### 'RandomEffectPredictors' — Names of columns in random-effects design matrix or cell arraycell array of length q | cell array of length R with elements of length q(r), r = 1, 2, ..., R

Names of columns in the random-effects design matrix or cell array Z, specified as the comma-separated pair consisting of 'RandomEffectPredictors' and either of the following:

• A cell array of length q when Z is an n-by-q design matrix. In this case, the default is {'z1','z2',...,'zQ'}.

• A cell array of length R, when Z is a cell array of length R with each element Z{r} of length q(r), r = 1, 2, ..., R. In this case, the default is {'z11','z12',...,'z1Q(1)'},...,{'zr1','zr2',...,'zrQ(r)'}.

For example, suppose you have correlated random effects for intercept and a variable named Acceleration. Then, you can specify the random-effects predictor names as follows.

Example: 'RandomEffectPredictors',{'Intercept','Acceleration'}

If you have two random effects terms, one for the intercept and the variable Acceleration grouped by variable g1, and the second for the intercept, grouped by the variable g2, then you specify the random-effects predictor names as follows.

Example: 'RandomEffectPredictors',{{'Intercept','Acceleration'},{'Intercept'}}

Data Types: cell

### 'ResponseVarName' — Name of response variable'y' (default) | string

Name of response variable, specified as the comma-separated pair consisting of 'ResponseVarName' and a string.

For example, if your response variable name is score, then you can specify it as follows.

Example: 'ResponseVarName','score'

Data Types: char

### 'RandomEffectGroups' — Names of random effects grouping variables'g' or {'g1','g2',...,'gR'} (default) | string | cell array of strings

Names of random effects grouping variables, specified as the comma-separated pair 'RandomEffectGroups' and either of the following:

• String — If there is only one random-effects term, that is, if G is a vector, then the value of 'RandomEffectGroups' is a string containing the name for the grouping variable G. The default is 'g'.

• Cell array of strings — If there are multiple random-effects terms, that is, if G is a cell array of length R, then the value of 'RandomEffectGroups' is a cell array of length R, where each cell contains the name for the grouping variable G{r}. The default is {'g1','g2',...,'gR'}.

For example, if you have two random-effects terms, z1 and z2, grouped by the grouping variables sex and subject, then you can specify the names of your grouping variables as follows.

Example: 'RandomEffectGroups',{'sex','subject'}

Data Types: char | cell

### 'CovariancePattern' — Pattern of covariance matrix'FullCholesky' (default) | string | square symmetric logical matrix | cell array of strings or logical matrices

Pattern of the covariance matrix of the random effects, specified as the comma-separated pair consisting of 'CovariancePattern' and a string, a square symmetric logical matrix, or a cell array of strings or logical matrices.

If there are R random-effects terms, then the value of 'CovariancePattern' must be a cell array of length R, where each element r of this cell array specifies the pattern of the covariance matrix of the random-effects vector associated with the rth random-effects term. The options for each element follow.

 'FullCholesky' Default. Full covariance matrix using the Cholesky parameterization. fitlme estimates all elements of the covariance matrix. 'Full' Full covariance matrix, using the log-Cholesky parameterization. fitlme estimates all elements of the covariance matrix. 'Diagonal' Diagonal covariance matrix. That is, off-diagonal elements of the covariance matrix are constrained to be 0. $\left(\begin{array}{ccc}{\sigma }_{b1}^{2}& 0& 0\\ 0& {\sigma }_{b2}^{2}& 0\\ 0& 0& {\sigma }_{b3}^{2}\end{array}\right)$ 'Isotropic' Diagonal covariance matrix with equal variances. That is, off-diagonal elements of the covariance matrix are constrained to be 0, and the diagonal elements are constrained to be equal. For example, if there are three random-effects terms with an isotropic covariance structure, this covariance matrix looks like$\left(\begin{array}{ccc}{\sigma }_{b}^{2}& 0& 0\\ 0& {\sigma }_{b}^{2}& 0\\ 0& 0& {\sigma }_{b}^{2}\end{array}\right)$where σ2b is the common variance of the random-effects terms. 'CompSymm' Compound symmetry structure. That is, common variance along diagonals and equal correlation between all random effects. For example, if there are three random-effects terms with a covariance matrix having a compound symmetry structure, this covariance matrix looks like$\left(\begin{array}{ccc}{\sigma }_{b1}^{2}& {\sigma }_{b1,b2}& {\sigma }_{b1,b2}\\ {\sigma }_{b1,b2}& {\sigma }_{b1}^{2}& {\sigma }_{b1,b2}\\ {\sigma }_{b1,b2}& {\sigma }_{b1,b2}& {\sigma }_{b1}^{2}\end{array}\right)$where σ2b1 is the common variance of the random-effects terms and σb1,b2 is the common covariance between any two random-effects term . PAT Square symmetric logical matrix. If 'CovariancePattern' is defined by the matrix PAT, and if PAT(a,b) = false, then the (a,b) element of the corresponding covariance matrix is constrained to be 0.

Example: 'CovariancePattern','Diagonal'

Example: 'CovariancePattern',{'Full','Diagonal'}

### 'FitMethod' — Method for estimating parameters'ML' (default) | 'REML'

Method for estimating parameters of the linear mixed-effects model, specified as the comma-separated pair consisting of 'FitMethod' and either of the following.

 'ML' Default. Maximum likelihood estimation 'REML' Restricted maximum likelihood estimation

Example: 'FitMethod','REML'

### 'Weights' — Observation weightsvector of scalar values

Observation weights, specified as the comma-separated pair consisting of 'Weights' and a vector of length n, where n is the number of observations.

Data Types: single | double

### 'Exclude' — Indices for rows to excludeuse all rows without NaNs (default) | vector of integer or logical values

Indices for rows to exclude from the linear mixed-effects model in the data, specified as the comma-separated pair consisting of 'Exclude' and a vector of integer or logical values.

For example, you can exclude the 13th and 67th rows from the fit as follows.

Example: 'Exclude',[13,67]

Data Types: single | double | logical

### 'DummyVarCoding' — Coding to use for dummy variables'reference' (default) | 'effects' | 'full'

Coding to use for dummy variables created from the categorical variables, specified as the comma-separated pair consisting of 'DummyVarCoding' and one of the following.

 'reference' Default. Coefficient for first category set to 0. 'effects' Coefficients sum to 0. 'full' One dummy variable for each category.

Example: 'DummyVarCoding','effects'

### 'Optimizer' — Optimization algorithm'quasinewton' (default) | 'fminunc'

Optimization algorithm, specified as the comma-separated pair consisting of 'Optimizer' and either of the following.

 'quasinewton' Default. Uses a trust region based quasi-Newton optimizer. Change the options of the algorithm using statset('LinearMixedModel'). If you don't specify the options, then LinearMixedModel uses the default options of statset('LinearMixedModel'). 'fminunc' You must have Optimization Toolbox™ to specify this option. Change the options of the algorithm using optimoptions('fminunc'). If you don't specify the options, then LinearMixedModel uses the default options of optimoptions('fminunc') with 'Algorithm' set to 'quasi-newton'.

Example: 'Optimizer','fminunc'

### 'OptimizerOptions' — Options for optimization algorithmstructure returned by statset | object returned by optimoptions

Options for the optimization algorithm, specified as the comma-separated pair consisting of 'OptimizerOptions' and a structure returned by statset('LinearMixedModel') or an object returned by optimoptions('fminunc').

• If 'Optimizer' is 'fminunc', then use optimoptions('fminunc') to change the options of the optimization algorithm. See optimoptions for the options 'fminunc' uses. If 'Optimizer' is 'fminunc' and you do not supply 'OptimizerOptions', then the default for LinearMixedModel is the default options created by optimoptions('fminunc') with 'Algorithm' set to 'quasi-newton'.

• If 'Optimizer' is 'quasinewton', then use statset('LinearMixedModel') to change the optimization parameters. If you don't change the optimization parameters, then LinearMixedModel uses the default options created by statset('LinearMixedModel'):

The 'quasinewton' optimizer uses the following fields in the structure created by statset('LinearMixedModel').

### 'TolFun' — Relative tolerance on gradient of objective function1e-6 (default) | positive scalar value

Relative tolerance on the gradient of the objective function, specified as a positive scalar value.

### 'TolX' — Absolute tolerance on step size1e-12 (default) | positive scalar value

Absolute tolerance on the step size, specified as a positive scalar value.

### 'MaxIter' — Maximum number of iterations allowed10000 (default) | positive scalar value

Maximum number of iterations allowed, specified as a positive scalar value.

### 'Display' — Level of display'off' (default) | 'iter' | 'final'

Level of display, specified as one of 'off', 'iter', or 'final'.

### 'StartMethod' — Method to start iterative optimization'default' (default) | 'random'

Method to start iterative optimization, specified as the comma-separated pair consisting of 'StartMethod' and either of the following.

 'default' Default. An internally defined default value. 'random' A random initial value.

Example: 'StartMethod','random'

### 'Verbose' — Indicator to display optimization process on screenfalse (default) | true

Indicator to display the optimization process on screen, specified as the comma-separated pair consisting of 'Verbose' and either false or true. Default is false.

The setting for 'Verbose' overrides the field 'Display' in 'OptimizerOptions'.

Example: 'Verbose',true

### 'CheckHessian' — Indicator to check positive definiteness of Hessianfalse (default) | true

Indicator to check the positive definiteness of the Hessian of the objective function with respect to unconstrained parameters at convergence, specified as the comma-separated pair consisting of 'CheckHessian' and either false or true. Default is false.

Specify 'CheckHessian' as true to verify optimality of the solution or to determine if the model is overparameterized in the number of covariance parameters.

Example: 'CheckHessian',true

## Output Arguments

expand all

### lme — Linear mixed-effects modelLinearMixedModel object

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

For properties and methods of this object, see LinearMixedModel.

## Alternative Functionality

You can also fit a linear mixed-effects model using fitlme(tbl,formula), where tbl is a table or dataset array containing the response y, the predictor variables X, and the grouping variables, and formula is of the form 'y ~ fixed + (random1|g1) + ... + (randomR|gR)'.

expand all

### Cholesky Parameterization

One of the assumptions of linear mixed-effects models is that the random effects have the following prior distribution.

$b~N\left(0,{\sigma }^{2}D\left(\theta \right)\right),$

where D is a q-by-q symmetric and positive semidefinite matrix, parameterized by a variance component vector θ, q is the number of variables in the random-effects term, and σ2 is the observation error variance. Since the covariance matrix of the random effects, D, is symmetric, it has q(q+1)/2 free parameters. Suppose L is the lower triangular Cholesky factor of D(θ) such that

$D\left(\theta \right)=L\left(\theta \right)L{\left(\theta \right)}^{T},$

then the q*(q+1)/2-by-1 unconstrained parameter vector θ is formed from elements in the lower triangular part of L.

For example, if

$L=\left[\begin{array}{ccc}{L}_{11}& 0& 0\\ {L}_{21}& {L}_{22}& 0\\ {L}_{31}& {L}_{32}& {L}_{33}\end{array}\right],$

then

$\theta =\left[\begin{array}{c}{L}_{11}\\ {L}_{21}\\ {L}_{31}\\ {L}_{22}\\ {L}_{32}\\ {L}_{33}\end{array}\right].$

### Log-Cholesky Parameterization

When the diagonal elements of L in Cholesky parameterization are constrained to be positive, then the solution for L is unique. Log-Cholesky parameterization is the same as Cholesky parameterization except that the logarithm of the diagonal elements of L are used to guarantee unique parameterization.

For example, for the 3-by-3 example in Cholesky parameterization, enforcing Lii ≥ 0,

$\theta =\left[\begin{array}{c}\mathrm{log}\left({L}_{11}\right)\\ {L}_{21}\\ {L}_{31}\\ \mathrm{log}\left({L}_{22}\right)\\ {L}_{32}\\ \mathrm{log}\left({L}_{33}\right)\end{array}\right].$