fitlmematrix

Fit linear mixed-effects model

Syntax

lme = fitlmematrix(X,y,Z,[])

lme = fitlmematrix(X,y,Z,G)

lme = fitlmematrix(___,Name,Value)

Description

lme= fitlmematrix(X,y,Z,[])creates a linear mixed-effects model of the responsesyusing the fixed-effects design matrixXand random-effects design matrix or matrices inZ.

[]implies that there is one group. That is, the grouping variableGisones(n,1), wherenis the number of observations. Usingfitlmematrix(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 designZand 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 responsesyusing the fixed-effects design matrixXand random-effects design matrixZor matrices inZ, and the grouping variable or variables inG.

example

lme= fitlmematrix(___,Name,Value)also creates a linear mixed-effects model with additional options specified by one or moreName,Valuepair 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

collapse all

No Grouping Variable Specified

Open Live Script

Load the sample data.

loadcarsmall

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 =毫升模型线性mixed-effects模型适合我nformation: 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 {'x1'} 43.575 2.3038 18.915 92 1.8371e-33 {'x2'} -0.0067097 0.0004242 -15.817 92 5.5373e-28 Lower Upper 39 48.151 -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 theZmatrix.

双([Z = Model_Year = = 70, Model_Year = = 76, Model_等Year==82]); lme = fitlmematrix(X,y,Z,[],'Covariancepattern','Isotropic')

lme =毫升模型线性mixed-effects模型适合我nformation: 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 {'x1'} 43.575 2.3038 18.915 92 1.8371e-33 {'x2'} -0.0067097 0.0004242 -15.817 92 5.5373e-28 Lower Upper 39 48.151 -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

Open Live Script

Load the sample data.

load('weight.mat');

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

DefineSubjectandProgramas 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_{i m} = β_{0} + β_{1} I W_{i} + β_{2} W e e k_{i} + β_{3} I {[P B]}_{i} + β_{4} I {[P C]}_{i} + β_{5} I {[P D]}_{i} \\ + β_{6} (W e e k_{i} * I {[P B]}_{i}) + β_{7} (W e e k_{i} * I {[P C]}_{i}) + β_{8} (W e e k_{i} * I {[P D]}_{i}) \\ + b_{0 m} + b_{1 m} W e e k_{i m} + ε_{i m}, \end{array}$

where $i$ = 1, 2, ..., 120, and $m$ = 1, 2, ..., 20. $β_{j}$ are the fixed-effects coefficients, $j$ = 0, 1, ..., 8, and $b_{0 m}$ and $b_{1 m}$ are random effects. $I W$ stands for initial weight and $I [\cdot]$ is a dummy variable representing a type of program. For example, $I [P B]_{i}$ is the dummy variable representing program type B. The random effects and observation error have the following prior distributions:

$b_{0 m} \sim N (0, σ_{0}^{2})$

$b_{1 m} \sim N (0, σ_{1}^{2})$

$ε_{i m} \sim N (0, σ^{2}) .$

Subject = nominal(Subject); Program = nominal(Program); D = dummyvar(Program);% Create dummy variables for ProgramX = [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 usingfitlmematrix与定义的德sign 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 =毫升模型线性mixed-effects模型适合我nformation: 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 {'Intercept' } 0.66105 0.25892 2.5531 111 0.012034 {'InitWeight'} 0.0031879 0.0013814 2.3078 111 0.022863 {'PrgB' } 0.36079 0.13139 2.746 111 0.0070394 {'PrgC' } -0.033263 0.13117 -0.25358 111 0.80029 {'PrgD' } 0.11317 0.13132 0.86175 111 0.39068 {'Week' } 0.1732 0.067454 2.5677 111 0.011567 {'Week_PrgB' } 0.038771 0.095394 0.40644 111 0.68521 {'Week_PrgC' } 0.030543 0.095394 0.32018 111 0.74944 {'Week_PrgD' } 0.033114 0.095394 0.34713 111 0.72915 Lower Upper 0.14798 1.1741 0.00045067 0.0059252 0.10044 0.62113 -0.29319 0.22666 -0.14706 0.3734 0.039536 0.30686 -0.15026 0.2278 -0.15849 0.21957 -0.15592 0.22214 Random effects covariance parameters (95% CIs): Group: Subject (20 Levels) Name1 Name2 Type Estimate {'Intercept'} {'Intercept'} {'std' } 0.18407 {'Week' } {'Intercept'} {'corr'} 0.66841 {'Week' } {'Week' } {'std' } 0.15033 Lower Upper 0.12281 0.27587 0.21076 0.88573 0.11004 0.20537 Group: Error Name Estimate Lower Upper {'Res Std'} 0.10261 0.087882 0.11981

Examine the fixed effects coefficients table. The row labeled'InitWeight'has a $p$ -value of 0.0228, and the row labeled'Week'has a $p$ -value of 0.0115. These $p$ -values 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 thecomparemethod.

Random Intercept Model

Open Live Script

Load the sample data.

loadflu

Thefludataset array has aDatevariable, 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).

适合一个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 estimateWtdILI, that shows which region each estimate is from, and the grouping variableDate.

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 byDate. The corresponding model is

$y_{i m} = β_{0} + β_{1} {W t d I L I}_{i m} + b_{0 m} + ε_{i m}, i = 1, 2, . . ., 468, m = 1, 2, . . ., 52,$

where $y_{i m}$ is the observation $i$ for level $m$ of grouping variableDate, $b_{0 m}$ is the random effect for level $m$ of the grouping variableDate, and $ε_{i m}$ is the observation error for observation $i$ . The random effect has the prior distribution,

$b_{0 m} \sim N (0, σ_{b}^{2}),$

and the error term has the distribution,

$ε_{i m} \sim 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 =毫升模型线性mixed-effects模型适合我nformation: 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 {'Intercept' } 0.16385 0.057525 2.8484 466 0.0045885 {'NationalRate'} 0.7236 0.032219 22.459 466 3.0502e-76 Lower Upper 0.050813 0.27689 0.66028 0.78691 Random effects covariance parameters (95% CIs): Group: Date (52 Levels) Name1 Name2 Type Estimate Lower {'Intercept'} {'Intercept'} {'std'} 0.17146 0.13227 Upper 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 $σ_{b}$ , 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 usingcomparemethod.

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} {y_{}^{ˆ}}_{28} = {β_{}^{ˆ}}_{0} + {β_{}^{ˆ}}_{1} {W t d I L I}_{28} + {b_{}^{ˆ}}_{10 / 30 / 2005} \\ = 0.1639 + 0.7236 * (1.343) + 0.3318 \\ = 1.46749, \end{array}$

where $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 STATSSTATS.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 thefitted(lme)method.

F = fitted(lme); F(28)

ans = 1.4674

Randomized Block Design

Open Live Script

Load the sample data.

load('shift.mat');

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,X1andX2, where

$Shift_Evening = {\begin{cases} 0, & if Morning \\ 1, & if Evening \\ - 1, & if Night \end{cases}$

$Shift_Morning = {\begin{cases} 1, & if Morning \\ 0, & if Evening \\ - 1, & if Night \end{cases}$

The model corresponds to

$\begin{array}{lrl} Morning Shift: & {QCDev}_{i m} & = β_{0} + β_{2} {Shift_Morning}_{i} + b_{0 m} + ϵ_{i m}, \\ Evening Shift: & {QCDev}_{i m} & = β_{0} + β_{1} {Shift_Evening}_{i} + b_{0 m} + ϵ_{i m}, \\ Night Shift: & {QCDev}_{i m} & = β_{0} - β_{1} {Shift_Evening}_{i} - β_{2} {Shift_Morning}_{i} + b_{0 m} + ϵ_{i m}, \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 the following distributions:

$b_{0 m} \sim N (0, σ_{b}^{2})$

and

$ε_{i m} \sim 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;

使用指定的符合线性mixed-effects模型ed 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 =线性mixed-effects模型适合REML模型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 {'Intercept'} 3.6525 0.94109 3.8812 12 0.0021832 {'S_Morning'} -0.91973 0.31206 -2.9473 12 0.012206 {'S_Evening'} -0.53293 0.31206 -1.7078 12 0.11339 Lower Upper 1.6021 5.703 -1.5997 -0.23981 -1.2129 0.14699 Random effects covariance parameters (95% CIs): Group: Operator (5 Levels) Name1 Name2 Type Estimate Lower {'Intercept'} {'Intercept'} {'std'} 2.0457 0.98207 Upper 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 =5×10.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} {y_{}^{ˆ}}_{Evening, Operator 3} = {β_{}^{ˆ}}_{0} + {β_{}^{ˆ}}_{1} Shift_Evening + {b_{}^{ˆ}}_{03} \\ = 3.6525 - 0.53293 - 2.1715 \\ = 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

Open Live Script

Load the sample data.

loadcarbig

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

${M P G}_{i m} = β_{0} + β_{1} {A c c}_{i} + β_{2} H P + b_{0 m} + {b_{1}}_{m} {A c c}_{i m} + ε_{i m}, m = 1, 2, 3,$

with the random-effects terms having the following prior distributions:

$b_{0 m} \sim N (0, σ_{0}^{2}),$

$b_{1 m} \sim N (0, σ_{1}^{2}),$

where $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 usingfitlmematrix与定义的德sign matrices and grouping variables.

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

lme =毫升模型线性mixed-effects模型适合我nformation: Number of observations 392 Fixed effects coefficients 3 Random effects coefficients 26 Covariance parameters 3 Formula: Linear Mixed Formula with 4 predictors. 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 {'Intercept' } 49.839 2.0518 24.291 389 {'Acceleration'} -0.58565 0.10846 -5.3995 389 {'Horsepower' } -0.16534 0.0071227 -23.213 389 pValue Lower Upper 5.6168e-80 45.806 53.873 1.1652e-07 -0.7989 -0.3724 1.9755e-75 -0.17934 -0.15133 Random effects covariance parameters (95% CIs): Group: Model_Year (13 Levels) Name1 Name2 Type Estimate Lower {'Intercept'} {'Intercept'} {'std'} 8.5771e-07 NaN Upper NaN Group: Model_Year (13 Levels) Name1 Name2 Type Estimate {'Acceleration'} {'Acceleration'} {'std'} 0.18783 Lower Upper 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} = (\begin{array}{l} b_{0 m} \\ b_{1 m} \end{array}) \sim N (0, (\begin{array}{cccccccccccccccccccc} σ_{0}^{2} & σ_{0, 1} \\ σ_{0, 1} & σ_{1}^{2} \end{array})),$

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 =毫升模型线性mixed-effects模型适合我nformation: Number of observations 392 Fixed effects coefficients 3 Random effects coefficients 26 Covariance parameters 4 Formula: Linear Mixed Formula with 4 predictors. 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 {'Intercept' } 50.133 2.2652 22.132 389 {'Acceleration'} -0.58327 0.13394 -4.3545 389 {'Horsepower' } -0.16954 0.0072609 -23.35 389 pValue Lower Upper 7.7727e-71 45.679 54.586 1.7075e-05 -0.84661 -0.31992 5.188e-76 -0.18382 -0.15527 Random effects covariance parameters (95% CIs): Group: Model_Year (13 Levels) Name1 Name2 Type Estimate {'Intercept' } {'Intercept' } {'std' } 3.3475 {'Acceleration'} {'Intercept' } {'corr'} -0.87971 {'Acceleration'} {'Acceleration'} {'std' } 0.33789 Lower Upper 1.2862 8.7119 -0.98501 -0.29675 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 thecomparemethod.

Specify the Covariance Pattern

Open Live Script

Load the sample data.

load('weight.mat');

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

DefineSubjectandProgramas 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_{i m} = β_{0} + β_{1} I W_{i} + β_{2} W e e k_{i} + β_{3} I {[P B]}_{i} + β_{4} I {[P C]}_{i} + β_{5} I {[P D]}_{i} \\ + b_{0 m} + b_{1 m} W e e k 2_{i m} + b_{2 m} W e e k 4_{i m} + b_{3 m} W e e k 6_{i m} + b_{4 m} W e e k 8_{i m} \\ + b_{5 m} W e e k 1 0_{i m} + b_{6 m} W e e k 1 2_{i m} + ε_{i m}, \end{array}$

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

$β_{j}$ are the fixed-effects coefficients, $j$ = 0, 1, ...,8, and $b_{0 m}$ and $b_{1 m}$ are random effects. $I W$ stands for initial weight and $I [.]$ is a dummy variable representing a type of program. For example, $I [P B]_{i}$ is the dummy variable representing program type B. The random effects and observation error have the following prior distributions:

$b_{0 m} \sim N (0, σ_{0}^{2})$

$b_{1 m} \sim N (0, σ_{1}^{2})$

$ε_{i m} \sim N (0, σ^{2}) .$

Fit the model usingfitlmematrix与定义的德sign 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 =毫升模型线性mixed-effects模型适合我nformation: Number of observations 120 Fixed effects coefficients 6 Random effects coefficients 40 Covariance parameters 2 Formula: Linear Mixed Formula with 7 predictors. 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 {'Intercept' } 0.4208 0.28169 1.4938 114 {'InitWeight'} 0.0045552 0.0015338 2.9699 114 {'PrgB' } 0.36993 0.12119 3.0525 114 {'PrgC' } -0.034009 0.1209 -0.28129 114 {'PrgD' } 0.121 0.12111 0.99911 114 {'Week' } 0.19881 0.037134 5.3538 114 pValue Lower Upper 0.13799 -0.13723 0.97883 0.0036324 0.0015168 0.0075935 0.0028242 0.12986 0.61 0.77899 -0.27351 0.2055 0.31986 -0.11891 0.36091 4.5191e-07 0.12525 0.27237 Random effects covariance parameters (95% CIs): Group: Subject (20 Levels) Name1 Name2 Type Estimate Lower {'Intercept'} {'Intercept'} {'std'} 0.16561 0.12896 Upper 0.21269 Group: Error Name Estimate Lower Upper {'Res Std'} 0.10272 0.088014 0.11987

Input Arguments

collapse all

`X`—Fixed-effects design matrix
n-by-pmatrix

Fixed-effects design matrix, specified as ann-by-pmatrix, wherenis the number of observations, andpis the number of fixed-effects predictor variables. Each row ofXcorresponds to one observation, and each column ofXcorresponds to one variable.

Data Types:single|double

`y`—Response values
n-by-1 vector

Response values, specified as ann-by-1 vector, wherenis the number of observations.

Data Types:single|double

`Z`—Random-effects design
n-by-qmatrix|cell array ofRn-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, thenZmust be ann-by-qmatrix, wherenis the number of observations andqis the number of variables in the random-effects term.
If there areRrandom-effects terms, thenZ必须是一个单元阵列的长度R. Each cell ofZcontains ann-by-q(r) design matrixZ{r},r= 1, 2, ...,R, corresponding to each random-effects term. Here,q(r) is the number of random effects term in therth random effects design matrix,Z{r}.

Data Types:single|double|cell

`G`—Grouping variable or variables
n-by-1 vector|cell array ofRn-by-1 vectors

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

If there is one random-effects term, thenGmust be ann-by-1 vector corresponding to a single grouping variable withMlevels or groups.
Gcan be a categorical vector, logical vector, numeric vector, character array, string array, or cell array of character vectors.
If there are multiple random-effects terms, thenG必须是一个单元阵列的长度R. Each cell ofGcontains a grouping variableG{r},r= 1, 2, ...,R, withM(r) levels.
G{r}can be a categorical vector, logical vector, numeric vector, character array, string array, or cell array of character vectors.

Name-Value Pair Arguments

Specify optional comma-separated pairs ofName,Valuearguments.Nameis the argument name andValueis the corresponding value.Namemust appear inside quotes. You can specify several name and value pair arguments in any order asName1,Value1,...,NameN,ValueN.

Example:'CovariancePattern','Diagonal','DummyVarCoding','full','Optimizer','fminunc'指定了一个随机协方差模式zero off-diagonal elements, creates a dummy variable for each level of a categorical variable, and uses thefminuncoptimization algorithm.

`'FixedEffectPredictors'`—Names of columns in fixed-effects design matrix
`{'x1','x2',...,'xP'}`(default) |string array or cell array of lengthp

Names of columns in the fixed-effects design matrixX, specified as the comma-separated pair consisting of'FixedEffectPredictors'and a string array or cell array of lengthp.

For example, if you have a constant term and two predictors, sayTimeSpentandGender, whereFemaleis the reference level forGender, as the fixed effects, then you can specify the names of your fixed effects in the following way.Gender_Malerepresents the dummy variable you must create for categoryMale. You can choose different names for these variables.

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

Data Types:string|cell

`'RandomEffectPredictors'`—Names of columns in random-effects design matrix or cell array
string array or cell array of lengthq|cell array of lengthRwith elements of lengthq(r),r= 1, 2, ...,R

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

A string array or cell array of lengthqwhenZis ann-by-qdesign matrix. In this case, the default is{'z1','z2',...,'zQ'}.
A cell array of lengthR, whenZis a cell array of lengthRwith each elementZ{r}of lengthq(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 namedAcceleration. 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 variableAccelerationgrouped by variableg1, and the second for the intercept, grouped by the variableg2, then you specify the random-effects predictor names as follows.

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

Data Types:string|cell

`'ResponseVarName'`—Name of response variable
`'y'`(default) |character vector|string scalar

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

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

Example:'ResponseVarName','score'

Data Types:char|string

`'RandomEffectGroups'`—Names of random effects grouping variables
`'g'`or`{'g1','g2',...,'gR'}`(default) |character vector|string scalar|string array|cell array of character vectors

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

Character vector or string scalar — If there is only one random-effects term, that is, ifGis a vector, then the value of'RandomEffectGroups'is the name for the grouping variableG. The default is'g'.
String array or cell array of character vectors — If there are multiple random-effects terms, that is, ifGis a cell array of lengthR, then the value of'RandomEffectGroups'is a string array or cell array of lengthR, where each element is the name for the grouping variableG{r}. The default is{'g1','g2',...,'gR'}.

For example, if you have two random-effects terms,z1andz2, grouped by the grouping variablessexandsubject, then you can specify the names of your grouping variables as follows.

Example:RandomEffectGroups,{“性”、“子ject'}

Data Types:char|string|cell

`'CovariancePattern'`—Pattern of covariance matrix
`'FullCholesky'`(default) |character vector|string scalar|square symmetric logical matrix|string array|cell array of character vectors or logical matrices

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

If there areRrandom-effects terms, then the value of'CovariancePattern'must be a string array or cell array of lengthR, where each elementrof the array specifies the pattern of the covariance matrix of the random-effects vector associated with therth 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. $(\begin{matrix} σ_{b 1}^{2} & 0 & 0 \\ 0 & σ_{b 2}^{2} & 0 \\ 0 & 0 & σ_{b 3}^{2} \end{matrix})$
`'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 $(\begin{matrix} σ_{b}^{2} & 0 & 0 \\ 0 & σ_{b}^{2} & 0 \\ 0 & 0 & σ_{b}^{2} \end{matrix})$ where σ²_bis 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 $(\begin{matrix} σ_{b 1}^{2} & σ_{b 1, b 2} & σ_{b 1, b 2} \\ σ_{b 1, b 2} & σ_{b 1}^{2} & σ_{b 1, b 2} \\ σ_{b 1, b 2} & σ_{b 1, b 2} & σ_{b 1}^{2} \end{matrix})$ where σ²_b1is the common variance of the random-effects terms and σ_b1,b2is 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'}

Data Types:char|string|logical|cell

`'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 weights
vector of scalar values

Observation weights, specified as the comma-separated pair consisting of'Weights'and a vector of lengthn, wherenis the number of observations.

Data Types:single|double

`'Exclude'`—Indices for rows to exclude
use 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 variables in this table.

Value	Description
`'reference'`(default)	`fitlmematrix`creates dummy variables with a reference group. This scheme treats the first category as a reference group and creates one less dummy variables than the number of categories. You can check the category order of a categorical variable by using the`categories`function, and change the order by using the`reordercats`function.
`'effects'`	`fitlmematrix`creates dummy variables using effects coding. This scheme uses –1 to represent the last category. This scheme creates one less dummy variables than the number of categories.
`'full'`	`fitlmematrix`creates full dummy variables. This scheme creates one dummy variable for each category.

For more details about creating dummy variables, seeAutomatic Creation of Dummy Variables.

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 algorithm
structure 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 bystatset('LinearMixedModel')or an object returned byoptimoptions('fminunc').

If'Optimizer'is'fminunc', then useoptimoptions('fminunc')to change the options of the optimization algorithm. Seeoptimoptionsfor the options'fminunc'uses. If'Optimizer'is'fminunc'and you do not supply'OptimizerOptions', then the default forLinearMixedModelis the default options created byoptimoptions('fminunc')with'Algorithm'set to'quasi-newton'.
If'Optimizer'is'quasinewton', then usestatset('LinearMixedModel')to change the optimization parameters. If you don’t change the optimization parameters, thenLinearMixedModeluses the default options created bystatset('LinearMixedModel'):

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

`TolFun`— Relative tolerance on gradient of objective function
`1e-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 size
`1e-12`(default) | positive scalar value

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

`MaxIter`— Maximum number of iterations allowed
`10000`(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.

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 eitherfalseortrue. Default isfalse.

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 eitherfalseortrue. Default isfalse.

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

Example:'CheckHessian',true

Output Arguments

collapse all

`lme`— Linear mixed-effects model
`LinearMixedModel`object

Linear mixed-effects model, returned as aLinearMixedModelobject.

More About

collapse all

Cholesky Parameterization

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

$b ~ N (0, σ^{2} D (θ)),$

whereDis aq-by-qsymmetric and positive semidefinite matrix, parameterized by a variance component vectorθ,qis the number of variables in the random-effects term, andσ²is the observation error variance. Since the covariance matrix of the random effects,D, is symmetric, it hasq(q+1)/2 free parameters. SupposeLis the lower triangular Cholesky factor ofD(θ) such that

$D (θ) = L (θ) L {(θ)}^{T},$

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

For example, if

$L = [\begin{matrix} L_{11} & 0 & 0 \\ L_{21} & L_{22} & 0 \\ L_{31} & L_{32} & L_{33} \end{matrix}],$

then

$θ = [\begin{matrix} L_{11} \\ L_{21} \\ L_{31} \\ L_{22} \\ L_{32} \\ L_{33} \end{matrix}] .$

Log-Cholesky Parameterization

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

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

$θ = [\begin{matrix} \log (L_{11}) \\ L_{21} \\ L_{31} \\ \log (L_{22}) \\ L_{32} \\ \log (L_{33}) \end{matrix}] .$

Alternative Functionality

You can also fit a linear mixed-effects model usingfitlme(tbl,formula), wheretblis a table or dataset array containing the responsey, the predictor variablesX, and the grouping variables, andformulais of the form'y ~ fixed + (random₁|g₁) + ... + (random_R|g_R)'.

fitlmematrix

Syntax

Description

Examples

No Grouping Variable Specified

Longitudinal Study with a Covariate

Random Intercept Model

Randomized Block Design

Correlated and Uncorrelated Random-Effects Terms

Specify the Covariance Pattern

Input Arguments

X—Fixed-effects design matrixn-by-pmatrix

y—Response valuesn-by-1 vector

Z—Random-effects designn-by-qmatrix|cell array ofRn-by-q(r) matrices,r= 1, 2, ...,R

G—Grouping variable or variablesn-by-1 vector|cell array ofRn-by-1 vectors

Name-Value Pair Arguments

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

'RandomEffectPredictors'—Names of columns in random-effects design matrix or cell arraystring array or cell array of lengthq|cell array of lengthRwith elements of lengthq(r),r= 1, 2, ...,R

'ResponseVarName'—Name of response variable'y'(default) |character vector|string scalar

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

'CovariancePattern'—Pattern of covariance matrix'FullCholesky'(default) |character vector|string scalar|square symmetric logical matrix|string array|cell array of character vectors or logical matrices

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

'Weights'—Observation weightsvector of scalar values

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

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

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

'OptimizerOptions'—Options for optimization algorithmstructure returned bystatset|object returned byoptimoptions

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

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

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

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

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

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

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

Output Arguments

lme— Linear mixed-effects modelLinearMixedModelobject

More About

Cholesky Parameterization

Log-Cholesky Parameterization

Alternative Functionality

See Also

Statistics and Machine Learning Toolbox Documentation

金宝app

Mastering Machine Learning: A Step-by-Step Guide with MATLAB

`X`—Fixed-effects design matrix
n-by-pmatrix

`y`—Response values
n-by-1 vector

`Z`—Random-effects design
n-by-qmatrix|cell array ofRn-by-q(r) matrices,r= 1, 2, ...,R

`G`—Grouping variable or variables
n-by-1 vector|cell array ofRn-by-1 vectors

`'FixedEffectPredictors'`—Names of columns in fixed-effects design matrix
`{'x1','x2',...,'xP'}`(default) |string array or cell array of lengthp

`'RandomEffectPredictors'`—Names of columns in random-effects design matrix or cell array
string array or cell array of lengthq|cell array of lengthRwith elements of lengthq(r),r= 1, 2, ...,R

`'ResponseVarName'`—Name of response variable
`'y'`(default) |character vector|string scalar

`'RandomEffectGroups'`—Names of random effects grouping variables
`'g'`or`{'g1','g2',...,'gR'}`(default) |character vector|string scalar|string array|cell array of character vectors

`'CovariancePattern'`—Pattern of covariance matrix
`'FullCholesky'`(default) |character vector|string scalar|square symmetric logical matrix|string array|cell array of character vectors or logical matrices

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

`'Weights'`—Observation weights
vector of scalar values

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

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

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

`'OptimizerOptions'`—Options for optimization algorithm
structure returned by`statset`|object returned by`optimoptions`

`TolFun`— Relative tolerance on gradient of objective function
`1e-6`(default) | positive scalar value

`TolX`— Absolute tolerance on step size
`1e-12`(default) | positive scalar value

`MaxIter`— Maximum number of iterations allowed
`10000`(default) | positive scalar value

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

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

`'Verbose'`—Indicator to display optimization process on screen
`false`(default) |`true`

`'CheckHessian'`—Indicator to check positive definiteness of Hessian
`false`(default) |`true`

`lme`— Linear mixed-effects model
`LinearMixedModel`object