Documentation Center

  • Trials
  • Product Updates

coefTest

Class: LinearMixedModel

Hypothesis test on fixed and random effects of linear mixed-effects model

Syntax

  • pVal = coefTest(lme) example
  • pVal = coefTest(lme,H) example
  • pVal = coefTest(lme,H,C) example
  • pVal = coefTest(lme,H,C,'REContrast',K) example
  • pVal = coefTest(___,Name,Value) example
  • [pVal,F,DF1,DF2] = coefTest(___) example

Description

example

pVal = coefTest(lme) returns the p-value for an F-test that all fixed-effects coefficients except for the intercept are 0.

example

pVal = coefTest(lme,H) returns the p-value for an F-test on fixed-effects coefficients of linear mixed-effects model lme, using the contrast matrix H. It tests the null hypothesis that H0: Hβ = 0, where β is the fixed-effects vector.

example

pVal = coefTest(lme,H,C) returns the p-value for an F-test on fixed-effects coefficients of the linear mixed-effects model lme, using the contrast matrix H. It tests the null hypothesis that H0: Hβ = C, where β is the fixed-effects vector.

example

pVal = coefTest(lme,H,C,'REContrast',K) returns the p-value for an F-test on the fixed- and random-effects coefficients of the linear mixed-effects model lme, using the contrast matrices H and K. It tests the null hypothesis that H0: Hβ + KB = C, where β is the fixed-effects vector and B is the random-effects vector.

example

pVal = coefTest(___,Name,Value) returns the p-value for an F-test on the fixed- and/or random-effects coefficients of the linear mixed-effects model lme, with additional options specified by one or more Name,Value pair arguments. It can accept any of the input arguments in the previous syntaxes.

For example, you can specify the method to compute the degrees of freedom.

example

[pVal,F,DF1,DF2] = coefTest(___) also returns the F-statistic F, and the numerator and denominator degrees of freedom for F, respectively DF1 and DF2.

Input Arguments

expand all

lme — Linear mixed-effects modelLinearMixedModel object

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

For properties and methods of this object, see LinearMixedModel.

H — Fixed-effects contrastsm-by-p matrix

Fixed-effects contrasts, specified as an m-by-p matrix, where p is the number of fixed-effects coefficients in lme. Each row of H represents one contrast. The columns of H (left to right) correspond to the rows of the p-by-1 fixed-effects vector beta (top to bottom), returned by the fixedEffects method.

Data Types: single | double

C — Hypothesized valuem-by-1 vector

Hypothesized value for testing the null hypothesis H*beta = C, specified as an m-by-1 matrix. Here, beta is the vector of fixed-effects estimates returned by the fixedEffects method.

Data Types: single | double

K — Random-effects contrastsm-by-q matrix

Random-effects contrasts, specified as an m-by-q matrix, where q is the number of random effects parameters in lme. The columns of K (left to right) correspond to the rows of the random-effects best linear unbiased predictor vector B (top to bottom), returned by the randomEffects method.

Data Types: single | double

Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside single quotes (' '). You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

'DFMethod' — Method for computing approximate denominator degrees of freedom'Residual' (default) | 'Satterthwaite' | 'None'

Method for computing the approximate denominator degrees of freedom for 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 np, 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

expand all

pValp-valuescalar value

p-value for the F-test on the fixed and/or random-effects coefficients of the linear mixed-effects model lme, returned as a scalar value.

FF-statisticscalar value

F-statistic, returned as a scalar value.

DF1 — Numerator degrees of freedom for Fscalar value

Numerator degrees of freedom for F, returned as a scalar value.

  • If you test the null hypothesis H0: Hβ = 0, or H0: Hβ = C, then DF1 is equal to the number of linearly independent rows in H.

  • If you test the null hypothesis H0: Hβ + KB= C, then DF1 is equal to the number of linearly independent rows in [H,K].

DF2 — Denominator degrees of freedom for Fscalar value

Denominator degrees of freedom for F, returned as a scalar value. The value of DF2 depends on the option you select for DFMethod.

Examples

expand all

Randomized Block Design

Navigate to a folder containing sample data.

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

Load the sample data.

load shift

The data shows the absolute deviations from the target quality characteristic measured from the products that five operators manufacture during three different 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 absolute 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 there is significant difference in the performance according to the time of the shift.

lme = fitlme(shift,'QCDev ~ Shift + (1|Operator)')
lme = 


Linear mixed-effects model fit by ML

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
    59.012    62.552    -24.506          49.012  

Fixed effects coefficients (95% CIs):
    Name                   Estimate    SE         tStat       DF    pValue       Lower      Upper  
    '(Intercept)'           3.1196     0.88681      3.5178    12    0.0042407     1.1874     5.0518
    'Shift_Morning'        -0.3868     0.48344    -0.80009    12      0.43921    -1.4401    0.66653
    'Shift_Night'           1.9856     0.48344      4.1072    12    0.0014535    0.93227     3.0389

Random effects covariance parameters (95% CIs):
Group: Operator (5 Levels)
    Name1                Name2                Type         Estimate    Lower      Upper 
    '(Intercept)'        '(Intercept)'        'std'        1.8297      0.94915    3.5272

Group: Error
    Name             Estimate    Lower      Upper 
    'Res Std'        0.76439     0.49315    1.1848

Test if all fixed-effects coefficients except for the intercept are 0.

pVal = coefTest(lme)
pVal =

   7.5956e-04

The small p-value indicates that not all fixed-effects coefficients are 0.

Test the significance of the Shift term using a contrast matrix.

H = [0 1 0; 0 0 1];
pVal = coefTest(lme,H)
pVal =

   7.5956e-04

Test the significance of the Shift term using the anova method.

anova(lme)
ans = 


    ANOVA marginal tests: DFMethod = 'Residual'

    Term                 FStat     DF1    DF2    pValue    
    '(Intercept)'        12.375    1      12      0.0042407
    'Shift'              13.864    2      12     0.00075956

The p-value for Shift, 0.00075956, is the same as the p-value of the previous hypothesis test.

Test if there is any difference between the evening and morning shifts.

pVal = coefTest(lme,[0 1 -1])
pVal =

   3.6147e-04

This small p-value indicates that the performance of the operators are not the same in the morning and the evening shifts.

Hypothesis Test for Fixed Effects

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 to 4 exercise programs, and their weight loss is recorded over six 2-week time periods. This is simulated data.

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

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

Fit a linear mixed-effects model where the initial weight, type of program, week, and the interaction between the week and type of program are the fixed effects. The intercept and week vary by subject.

lme = fitlme(tbl,'y ~ InitialWeight + Program*Week + (Week|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:
    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.66105      0.25892      2.5531    111     0.012034       0.14798       1.1741
    'InitialWeight'         0.0031879    0.0013814      2.3078    111     0.022863    0.00045067    0.0059252
    'Program_B'               0.36079      0.13139       2.746    111    0.0070394       0.10044      0.62113
    'Program_C'             -0.033263      0.13117    -0.25358    111      0.80029      -0.29319      0.22666
    'Program_D'               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
    'Program_B:Week'         0.038771     0.095394     0.40644    111      0.68521      -0.15026       0.2278
    'Program_C:Week'         0.030543     0.095394     0.32018    111      0.74944      -0.15849      0.21957
    'Program_D:Week'         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

Test for the significance of the interaction between Program and Week.

H = [0 0 0 0 0 0 1 0 0; 
     0 0 0 0 0 0 0 1 0;
     0 0 0 0 0 0 0 0 1];
pVal = coefTest(lme,H)
pVal =

    0.9775

The high p-value indicates that the interaction between Program and Week is not statistically significant.

Now, test whether all coefficients involving Program are 0.

H = [0 0 1 0 0 0 0 0 0;
     0 0 0 1 0 0 0 0 0;
     0 0 0 0 1 0 0 0 0;
     0 0 0 0 0 0 1 0 0; 
     0 0 0 0 0 0 0 1 0;
     0 0 0 0 0 0 0 0 1];
C = [0;0;0;0;0;0];
pVal = coefTest(lme,H,C)
pVal =

    0.0274

The p-value of 0.0274 indicates that not all coefficients involving Program are zero.

Test for Fixed and Random Effects

Load the sample data.

load flu

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

To fit a linear-mixed effects model, your data must be in a properly formatted dataset array. To fit a linear mixed-effects model with the influenza rates as the responses and region as the predictor variable, combine the nine columns corresponding to the regions into a tall array. The new dataset array, flu2, must have the response variable, FluRate, the nominal variable, Region, 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);

Fit a linear mixed-effects model with fixed effects for the region and a random intercept that varies by Date.

lme = fitlme(flu2,'FluRate ~ 1 + Region + (1|Date)')
Linear mixed-effects model fit by ML

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

Formula:
    FluRate ~ 1 + Region + (1|Date)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    318.71    364.35    -148.36          296.71  

Fixed effects coefficients (95% CIs):
    Name                      Estimate    SE          tStat      DF     pValue        Lower        Upper    
    '(Intercept)'               1.2233    0.096678     12.654    459     1.085e-31       1.0334       1.4133
    'Region_MidAtl'           0.010192    0.052221    0.19518    459       0.84534    -0.092429      0.11281
    'Region_ENCentral'        0.051923    0.052221     0.9943    459        0.3206    -0.050698      0.15454
    'Region_WNCentral'         0.23687    0.052221     4.5359    459    7.3324e-06      0.13424      0.33949
    'Region_SAtl'             0.075481    0.052221     1.4454    459       0.14902     -0.02714       0.1781
    'Region_ESCentral'         0.33917    0.052221      6.495    459    2.1623e-10      0.23655      0.44179
    'Region_WSCentral'           0.069    0.052221     1.3213    459       0.18705    -0.033621      0.17162
    'Region_Mtn'              0.046673    0.052221    0.89377    459       0.37191    -0.055948      0.14929
    'Region_Pac'              -0.16013    0.052221    -3.0665    459     0.0022936     -0.26276    -0.057514

Random effects covariance parameters (95% CIs):
Group: Date (52 Levels)
    Name1                Name2                Type         Estimate    Lower     Upper  
    '(Intercept)'        '(Intercept)'        'std'        0.6443      0.5297    0.78368

Group: Error
    Name             Estimate    Lower      Upper
    'Res Std'        0.26627     0.24878    0.285

Test the hypothesis that the random effects-term for week 10/9/2005 is zero.

[~,~,STATS] = randomEffects(lme); % Compute the random-effects statistics (STATS)
STATS.Level = nominal(STATS.Level);
K = zeros(length(STATS),1);
K(STATS.Level == '10/9/2005') = 1;
pVal = coefTest(lme,[0 0 0 0 0 0 0 0 0],0,'REContrast',K')
pVal =

    0.1692

Refit the model this time with a random intercept and slope.

lme = fitlme(flu2,'FluRate ~ 1 + Region + (1 + Region|Date)');

Test the hypothesis that the combined coefficient of region WNCentral for week 10/9/2005 is zero.

[~,~,STATS] = randomEffects(lme); STATS.Level = nominal(STATS.Level);
K = zeros(length(STATS),1);
K(STATS.Level == '10/9/2005' & flu2.Region == 'WNCentral') = 1;
pVal = coefTest(lme,[0 0 0 1 0 0 0 0 0],0,'REContrast',K')
ans =

   1.4536e-12

Also return the F-statistic with the numerator and denominator degrees of freedom.

[pVal,F,DF1,DF2] = coefTest(lme,[0 0 0 1 0 0 0 0 0],0,'REContrast',K')
pVal =

   1.4536e-12


F =

   53.0081


DF1 =

     1


DF2 =

   459

Repeat the test using the Satterthwaite approximation for the denominator degrees of freedom.

[pVal,F,DF1,DF2] = coefTest(lme,[0 0 0 1 0 0 0 0 0],0,'REContrast',K','DFMethod','Satterthwaite')
pVal =

    0.0055


F =

   53.0081


DF1 =

     1


DF2 =

    2.9795

See Also

| |

Was this topic helpful?