anova

Analysis of variance for linear mixed-effects model

Syntax

stats = anova(lme)

stats = anova(lme,Name,Value)

Description

stats = anova(lme) returns the dataset array stats that includes the results of the F-tests for each fixed-effects term in the linear mixed-effects model lme.

stats = anova(lme,Name,Value) also returns the dataset array stats with additional options specified by one or more Name,Value pair arguments.

Examples

F-Tests for Fixed Effects

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 deviation of the quality characteristics from the target value. This is simulated data.

Shift and Operator are nominal variables.

shift.Shift = nominal(shift.Shift);
shift.Operator = nominal(shift.Operator);

Fit a linear mixed-effects model with a random intercept grouped by operator to assess if performance significantly differs according to the time of the shift. Use the restricted maximum likelihood method and 'effects' contrasts.

'effects' contrasts indicate that the coefficients sum to 0, and fitlme creates two contrast-coded variables in the fixed-effects design matrix, $X$1 and $X$2, where

$$Shift\_Evening=\{\begin{array}{c}0,ifMorning\\ 1,ifEvening\\ -1,ifNight\end{array}$$

and

$$Shift\_Morning=\{\begin{array}{c}1,ifMorning\\ 0,ifEvening\\ -1,ifNight\end{array}.$$

The model corresponds to

$$\begin{array}{l}MorningShift:QCDe{v}_{im}={\beta }_0+{\beta }_{2}Shift\_Mornin{g}_i+b_{0m}+{\epsilon }_{im},m=1,2,...,5,\\ EveningShift:QCDe{v}_{im}={\beta }_0+{\beta }_{1}Shift\_Evenin{g}_i+b_{0m}+{\epsilon }_{im},\\ NightShift:QCDe{v}_{im}={\beta }_0-{\beta }_{1}Shift\_Evenin{g}_i-{\beta }_{2}Shift\_Mornin{g}_i+b_{0m}+{\epsilon }_{im},\end{array}$$

where $$b$$ ~ N(0, $${\sigma }_b^2$$ ) and $$ϵ$$ ~ N(0, $${\sigma }^2$$ ).

lme = fitlme(shift,'QCDev ~ Shift + (1|Operator)',...
'FitMethod','REML','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:
    QCDev ~ 1 + Shift + (1 | 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
    {'Shift_Evening'}        -0.53293    0.31206    -1.7078    12      0.11339    -1.2129     0.14699
    {'Shift_Morning'}        -0.91973    0.31206    -2.9473    12     0.012206    -1.5997    -0.23981

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

Perform an $$F$$-test to determine if all fixed-effects coefficients are 0.

anova(lme)

ans = 
    ANOVA MARGINAL TESTS: DFMETHOD = 'RESIDUAL'

    Term                   FStat     DF1    DF2    pValue   
    {'(Intercept)'}        15.063    1      12     0.0021832
    {'Shift'      }        11.091    2      12     0.0018721

The $$p$$-value for the constant term, 0.0021832, is the same as in the coefficient table in the lme display. The $$p$$-value of 0.0018721 for Shift measures the combined significance for both coefficients representing Shift.

ANOVA for Fixed-Effects in LME Model

Load the sample data.

load('fertilizer.mat')

The dataset array includes data from a split-plot experiment, where soil is divided into three blocks based on the soil type: sandy, silty, and loamy. Each block is divided into five plots, where five types of tomato plants (cherry, heirloom, grape, vine, and plum) are randomly assigned to these plots. The tomato plants in the plots are then divided into subplots, where each subplot is treated by one of four fertilizers. This is simulated data.

Store the data in a dataset array called ds, for practical purposes, and define Tomato, Soil, and Fertilizer as categorical variables.

ds = fertilizer;
ds.Tomato = nominal(ds.Tomato);
ds.Soil = nominal(ds.Soil);
ds.Fertilizer = nominal(ds.Fertilizer);

Fit a linear mixed-effects model, where Fertilizer and Tomato are the fixed-effects variables, and the mean yield varies by the block (soil type) and the plots within blocks (tomato types within soil types) independently. Use the 'effects' contrasts when fitting the data for the type III sum of squares.

lme = fitlme(ds,'Yield ~ Fertilizer * Tomato + (1|Soil) + (1|Soil:Tomato)',...
'DummyVarCoding','effects')

lme = 
Linear mixed-effects model fit by ML

Model information:
    Number of observations              60
    Fixed effects coefficients          20
    Random effects coefficients         18
    Covariance parameters                3

Formula:
    Yield ~ 1 + Tomato*Fertilizer + (1 | Soil) + (1 | Soil:Tomato)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    522.57    570.74    -238.29          476.57  

Fixed effects coefficients (95% CIs):
    Name                                    Estimate    SE        tStat       DF    pValue        Lower      Upper  
    {'(Intercept)'                 }           104.6    3.3008       31.69    40    5.9086e-30     97.929     111.27
    {'Tomato_Cherry'               }             1.4    5.9353     0.23588    40       0.81473    -10.596     13.396
    {'Tomato_Grape'                }         -7.7667    5.9353     -1.3085    40       0.19816    -19.762     4.2291
    {'Tomato_Heirloom'             }         -11.183    5.9353     -1.8842    40      0.066821    -23.179    0.81242
    {'Tomato_Plum'                 }          30.233    5.9353      5.0938    40     8.777e-06     18.238     42.229
    {'Fertilizer_1'                }         -28.267    2.3475     -12.041    40    7.0265e-15    -33.011    -23.522
    {'Fertilizer_2'                }         -1.9333    2.3475    -0.82356    40       0.41507    -6.6779     2.8112
    {'Fertilizer_3'                }          10.733    2.3475      4.5722    40     4.577e-05     5.9888     15.478
    {'Tomato_Cherry:Fertilizer_1'  }        -0.73333    4.6951    -0.15619    40       0.87667    -10.222     8.7558
    {'Tomato_Grape:Fertilizer_1'   }         -7.5667    4.6951     -1.6116    40       0.11491    -17.056     1.9224
    {'Tomato_Heirloom:Fertilizer_1'}          5.1833    4.6951       1.104    40       0.27619    -4.3058     14.672
    {'Tomato_Plum:Fertilizer_1'    }          2.7667    4.6951     0.58927    40       0.55899    -6.7224     12.256
    {'Tomato_Cherry:Fertilizer_2'  }             7.6    4.6951      1.6187    40       0.11337    -1.8891     17.089
    {'Tomato_Grape:Fertilizer_2'   }            -1.9    4.6951    -0.40468    40       0.68787    -11.389     7.5891
    {'Tomato_Heirloom:Fertilizer_2'}          5.5167    4.6951       1.175    40       0.24695    -3.9724     15.006
    {'Tomato_Plum:Fertilizer_2'    }            -3.9    4.6951    -0.83066    40        0.4111    -13.389     5.5891
    {'Tomato_Cherry:Fertilizer_3'  }         -6.0667    4.6951     -1.2921    40       0.20373    -15.556     3.4224
    {'Tomato_Grape:Fertilizer_3'   }          3.7667    4.6951     0.80226    40       0.42714    -5.7224     13.256
    {'Tomato_Heirloom:Fertilizer_3'}          3.1833    4.6951     0.67802    40       0.50167    -6.3058     12.672
    {'Tomato_Plum:Fertilizer_3'    }             1.1    4.6951     0.23429    40       0.81596    -8.3891     10.589

Random effects covariance parameters (95% CIs):
Group: Soil (3 Levels)
    Name1                  Name2                  Type           Estimate    Lower       Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        2.5028      0.027711    226.05

Group: Soil:Tomato (15 Levels)
    Name1                  Name2                  Type           Estimate    Lower     Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        10.225      6.1497    17.001

Group: Error
    Name               Estimate    Lower     Upper 
    {'Res Std'}        10.499      8.5389    12.908

Perform an analysis of variance to test for the fixed-effects.

anova(lme)

ans = 
    ANOVA MARGINAL TESTS: DFMETHOD = 'RESIDUAL'

    Term                         FStat     DF1    DF2    pValue    
    {'(Intercept)'      }        1004.2     1     40     5.9086e-30
    {'Tomato'           }        7.1663     4     40     0.00018935
    {'Fertilizer'       }        58.833     3     40     1.0024e-14
    {'Tomato:Fertilizer'}        1.4182    12     40        0.19804

The $$p$$-value for the constant term, 5.9086e-30, is the same as in the coefficient table in the lme display. The $$p$$-values of 0.00018935, 1.0024e-14, and 0.19804 for Tomato, Fertilizer, and Tomato:Fertilizer represent the combined significance for all tomato coefficients, fertilizer coefficients, and coefficients representing the interaction between the tomato and fertilizer, respectively. The $$p$$-value of 0.19804 indicates that the interaction between tomato and fertilizer is not significant.

Satterthwaite Approximation for Degrees of Freedom

Load the sample data.

load('weight.mat')

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.

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

tbl = table(InitialWeight,Program,Subject,Week,y);
tbl.Subject = nominal(tbl.Subject);
tbl.Program = nominal(tbl.Program);

Fit the model using the 'effects' contrasts.

lme = fitlme(tbl,'y ~ InitialWeight + Program*Week + (Week|Subject)',...
		'DummyVarCoding','effects')

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:
    y ~ 1 + InitialWeight + Program*Week + (1 + Week | Subject)

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.77122      0.24309      3.1725    111     0.0019549       0.28951       1.2529
    {'InitialWeight' }        0.0031879    0.0013814      2.3078    111      0.022863    0.00045067    0.0059252
    {'Program_A'     }         -0.11017     0.080377     -1.3707    111       0.17323      -0.26945       0.0491
    {'Program_B'     }          0.25061      0.08045      3.1151    111     0.0023402      0.091195      0.41003
    {'Program_C'     }         -0.14344     0.080475     -1.7824    111      0.077424       -0.3029     0.016031
    {'Week'          }          0.19881     0.033727      5.8946    111    4.1099e-08       0.13198      0.26564
    {'Program_A:Week'}        -0.025607     0.058417    -0.43835    111       0.66198      -0.14136     0.090149
    {'Program_B:Week'}         0.013164     0.058417     0.22535    111       0.82212      -0.10259      0.12892
    {'Program_C:Week'}        0.0049357     0.058417    0.084492    111       0.93282      -0.11082      0.12069

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.21077    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.022863 and 4.1099e-08 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 that are in program A. The lower and upper limits of the covariance parameters for the random effects do not include zero, thus they are significant.

Perform an F-test that all fixed-effects coefficients are zero.

anova(lme)

ans = 
    ANOVA MARGINAL TESTS: DFMETHOD = 'RESIDUAL'

    Term                     FStat       DF1    DF2    pValue    
    {'(Intercept)'  }          10.065    1      111     0.0019549
    {'InitialWeight'}           5.326    1      111      0.022863
    {'Program'      }          3.6798    3      111      0.014286
    {'Week'         }          34.747    1      111    4.1099e-08
    {'Program:Week' }        0.066648    3      111       0.97748

The $$p$$-values for the constant term, initial weight, and week are the same as in the coefficient table in the previous lme output display. The $$p$$-value of 0.014286 for Program represents the combined significance for all program coefficients. Similarly, the $$p$$-value for the interaction between program and week (Program:Week) measures the combined significance for all coefficients representing this interaction.

Now, use the Satterthwaite method to compute the degrees of freedom.

anova(lme,'DFMethod','satterthwaite')

ans = 
    ANOVA MARGINAL TESTS: DFMETHOD = 'SATTERTHWAITE'

    Term                     FStat       DF1    DF2       pValue    
    {'(Intercept)'  }          10.065    1      20.445      0.004695
    {'InitialWeight'}           5.326    1          20      0.031827
    {'Program'      }          3.6798    3       19.14      0.030233
    {'Week'         }          34.747    1          20    9.1346e-06
    {'Program:Week' }        0.066648    3          20       0.97697

The Satterthwaite method produces smaller denominator degrees of freedom and slightly larger $$p$$-values.

Input Arguments

lme — Linear mixed-effects model

LinearMixedModel object

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

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Example: stats = anova(lme,'DFMethod','satterthwaite')

DFMethod — Method for computing approximate degrees of freedom

'residual' (default) | 'satterthwaite' | 'none'

Method for computing approximate degrees of freedom to use in the F-test, specified as the comma-separated pair consisting of 'DFMethod' and one of the following.

'residual'	Default. The degrees of freedom are assumed to be constant and equal to n – p, where n is the number of observations and p is the number of fixed effects.
'satterthwaite'	Satterthwaite approximation.
'none'	All degrees of freedom are set to infinity.

For example, you can specify the Satterthwaite approximation as follows.

Example: 'DFMethod','satterthwaite'

Output Arguments

stats — Results of F-tests for fixed-effects terms

dataset array

Results of F-tests for fixed-effects terms, returned as a dataset array with the following columns.

Term	Name of the fixed effects term
Fstat	F-statistic for the term
DF1	Numerator degrees of freedom for the F-statistic
DF2	Denominator degrees of freedom for the F-statistic
pValue	p-value of the test for the term

There is one row for each fixed-effects term. Each term is a continuous variable, a grouping variable, or an interaction between two or more continuous or grouping variables. For each fixed-effects term, anova performs an F-test (marginal test) to determine if all coefficients representing the fixed-effects term are 0. To perform tests for the type III hypothesis, you must use the 'effects' contrasts while fitting the linear mixed-effects model.

Tips

• For each fixed-effects term, anova performs an F-test (marginal test), that all coefficients representing the fixed-effects term are 0. To perform tests for type III hypotheses, you must set the 'DummyVarCoding' name-value pair argument to 'effects' contrasts while fitting your linear mixed-effects model.

Version History

Introduced in R2013b