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.

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

**MathWorks Machine Translation**

The automated translation of this page is provided by a general purpose third party translator tool.

MathWorks does not warrant, and disclaims all liability for, the accuracy, suitability, or fitness for purpose of the translation.

Fit linear mixed-effects model

creates
a linear mixed-effects model of the responses `lme`

= fitlmematrix(`X`

,`y`

,`Z`

,[])`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.

also
creates a linear mixed-effects model with additional options specified
by one or more `lme`

= fitlmematrix(___,`Name,Value`

)`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.

Load the sample data.

`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

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 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{\hspace{1em}}\text{\hspace{1em}}+{\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{\hspace{1em}}\text{\hspace{1em}}+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 *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 type B. The random effects
and observation error have these prior distributions: *b*_{0m}~N(0,σ^{2}_{0}), *b*_{1m}~N(0,σ^{2}_{1}),
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.

Load the sample data.

`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 an 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{\hspace{1em}}i=1,2,\mathrm{...},468,\text{\hspace{1em}}m=1,2,\mathrm{...},52,$$

where *y*_{im} is
the observation *i* for level *m* of
grouping variable `Date`

. *b*_{0m} 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, *b*_{0m} ~
N(0,σ^{2}_{FR})
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 σ^{2}_{b0m},
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}{\widehat{y}}_{28}={\widehat{\beta}}_{0}+{\widehat{\beta}}_{1}WtdIL{I}_{28}+{\widehat{b}}_{10/30/2005}\\ \text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}=0.1639+0.7236*(1.343)+0.3318\\ \text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}=1.46749,\end{array}$$

where $$\widehat{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

Navigate to a folder containing sample data.

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

Load the sample data.

`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, *X*1
and *X*2, where

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

The model corresponds to

$$\begin{array}{l}\text{MorningShift:}QCDe{v}_{im}={\beta}_{0}+{\beta}_{2}Shift\text{\_}Mornin{g}_{i}+{b}_{0m}+{\epsilon}_{im},\text{\hspace{1em}}m=1,2,\mathrm{...},5,\\ \text{EveningShift:}QCDe{v}_{im}={\beta}_{0}+{\beta}_{1}Shift\text{\_}Evenin{g}_{i}+{b}_{0m}+{\epsilon}_{im},\\ \text{NightShift:}QCDe{v}_{im}={\beta}_{0}-{\beta}_{1}Shift\text{\_}Evenin{g}_{i}-{\beta}_{2}Shift\text{\_}Mornin{g}_{i}+{b}_{0m}+{\epsilon}_{im},\end{array}$$

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: *b*_{0m} ~
N(0,σ^{2}_{b0m})
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}{\widehat{y}}_{Evening,Operator3}={\widehat{\beta}}_{0}+{\widehat{\beta}}_{1}Shift\text{\_}Evening+{\widehat{b}}_{03}\\ \text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{0.17em}}=3.6525-0.53293-2.1715\\ \text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{0.17em}}=\mathrm{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

Load the sample data.

`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{\hspace{1em}}m=1,2,3,$$

with the random-effects
terms having these distributions: *b*_{0m} ~
N(0,σ^{2}_{0}),
and *b*_{1m} ~
N(0,σ^{2}_{1}). *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.

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 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{\hspace{1em}}\text{\hspace{1em}}+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{\hspace{1em}}\text{\hspace{1em}}+{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 *b*_{1m} 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 type B. The random effects
and observation error have these prior distributions: *b*_{0m}~N(0,σ^{2}_{0}), *b*_{1m}~N(0,σ^{2}_{1}),
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

`X`

— Fixed-effects design matrixFixed-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 valuesResponse values, specified as an *n*-by-1 vector,
where *n* is the number of observations.

**Data Types: **`single`

| `double`

`Z`

— Random-effects designRandom-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*r*th random effects design matrix,`Z{r}`

.

**Data Types: **`single`

| `double`

| `cell`

`G`

— Grouping variable or variablesGrouping 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 character vectors.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 character vectors.

**Data Types: **`single`

| `double`

| `char`

| `cell`

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`

.

`'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 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

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) | character vectorName of response variable, specified as the comma-separated
pair consisting of `'ResponseVarName'`

and a character
vector.

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) | character vector | cell array of character vectorsNames of random effects grouping variables, specified as the
comma-separated pair `'RandomEffectGroups'`

and either
of the following:

Character vector — If there is only one random-effects term, that is, if

`G`

is a vector, then the value of`'RandomEffectGroups'`

is a character vector containing the name for the grouping variable`G`

. The default is`'g'`

.Cell array of character vectors — 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) | character vector | square symmetric logical matrix | cell array of character vectors or logical matricesPattern of the covariance matrix of the random effects, specified
as the comma-separated pair consisting of `'CovariancePattern'`

and
a character vector, a square symmetric logical matrix, or a cell array
of character vectors 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 *r*th
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 σ |

`'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)$$ ^{2}_{b1} 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 valuesIndices 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.

Value | Description |
---|---|

`'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 function`1e-6`

(default) | positive scalar valueRelative tolerance on the gradient of the objective function, specified as a positive scalar value.

`'TolX'`

— Absolute tolerance on step size`1e-12`

(default) | positive scalar valueAbsolute tolerance on the step size, specified as a positive scalar value.

`'MaxIter'`

— Maximum number of iterations allowed`10000`

(default) | positive scalar valueMaximum 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.

Value | Description |
---|---|

`'default'` | An internally defined default value |

`'random'` | A random initial value |

**Example: **`'StartMethod','random'`

`'Verbose'`

— Indicator to display optimization process on screen`false`

(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 Hessian`false`

(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`

`lme`

— Linear mixed-effects model`LinearMixedModel`

objectLinear mixed-effects model, returned as a `LinearMixedModel`

object.

For properties and methods of this object, see `LinearMixedModel`

.

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 + (random`

. _{1}|g_{1})
+ ... + (random_{R}|g_{R})'

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].$$

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 *L*_{ii} ≥
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].$$

Was this topic helpful?

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.

You can also select a location from the following list:

- Canada (English)
- United States (English)

- 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)