coefTest

Class: LinearMixedModel

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

expand all in page

Syntax

pVal = coefTest(lme)
pVal = coefTest(lme,H)
pVal = coefTest(lme,H,C)
pVal = coefTest(lme,H,C,Name,Value)
[pVal,F,DF1,DF2] = coefTest(___)

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,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. For example, 'REContrast',K tells coefTest to test the null hypothesis that H0: Hβ + KB = C, where β is the fixed-effects vector and B is the random-effects vector.

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

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

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

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

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 quotes. You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

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

Random-effects contrasts, specified as the comma-separated pair consisting of 'REContrast' and an m-by-q matrix K, 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

Output Arguments

expand all

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.

F-statistic, returned as a scalar 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].

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

Open Live Script

Load the sample data. 

load('shift.mat')

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   
    {'(Intercept)'  }         3.1196     0.88681      3.5178    12    0.0042407
    {'Shift_Morning'}        -0.3868     0.48344    -0.80009    12      0.43921
    {'Shift_Night'  }         1.9856     0.48344      4.1072    12    0.0014535


    Lower      Upper  
     1.1874     5.0518
    -1.4401    0.66653
    0.93227     3.0389

Random effects covariance parameters (95% CIs):
Group: Operator (5 Levels)
    Name1                  Name2                  Type           Estimate
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        1.8297  


    Lower      Upper 
    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.

Open Live Script

Load the sample data. 

load('weight.mat')

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 
    {'(Intercept)'   }          0.66105      0.25892      2.5531    111
    {'InitialWeight' }        0.0031879    0.0013814      2.3078    111
    {'Program_B'     }          0.36079      0.13139       2.746    111
    {'Program_C'     }        -0.033263      0.13117    -0.25358    111
    {'Program_D'     }          0.11317      0.13132     0.86175    111
    {'Week'          }           0.1732     0.067454      2.5677    111
    {'Program_B:Week'}         0.038771     0.095394     0.40644    111
    {'Program_C:Week'}         0.030543     0.095394     0.32018    111
    {'Program_D:Week'}         0.033114     0.095394     0.34713    111


    pValue       Lower         Upper    
     0.012034       0.14798       1.1741
     0.022863    0.00045067    0.0059252
    0.0070394       0.10044      0.62113
      0.80029      -0.29319      0.22666
      0.39068      -0.14706       0.3734
     0.011567      0.039536      0.30686
      0.68521      -0.15026       0.2278
      0.74944      -0.15849      0.21957
      0.72915      -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

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.

Open Live Script

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 an 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)')
lme = 
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 
    {'(Intercept)'     }          1.2233    0.096678     12.654    459
    {'Region_MidAtl'   }        0.010192    0.052221    0.19518    459
    {'Region_ENCentral'}        0.051923    0.052221     0.9943    459
    {'Region_WNCentral'}         0.23687    0.052221     4.5359    459
    {'Region_SAtl'     }        0.075481    0.052221     1.4454    459
    {'Region_ESCentral'}         0.33917    0.052221      6.495    459
    {'Region_WSCentral'}           0.069    0.052221     1.3213    459
    {'Region_Mtn'      }        0.046673    0.052221    0.89377    459
    {'Region_Pac'      }        -0.16013    0.052221    -3.0665    459


    pValue        Lower        Upper    
     1.085e-31       1.0334       1.4133
       0.84534    -0.092429      0.11281
        0.3206    -0.050698      0.15454
    7.3324e-06      0.13424      0.33949
       0.14902     -0.02714       0.1781
    2.1623e-10      0.23655      0.44179
       0.18705    -0.033621      0.17162
       0.37191    -0.055948      0.14929
     0.0022936     -0.26276    -0.057514

Random effects covariance parameters (95% CIs):
Group: Date (52 Levels)
    Name1                  Name2                  Type           Estimate
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        0.6443  


    Lower     Upper  
    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')
pVal = 1.1939e-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.1939e-12

F = 53.4396

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 = NaN

F = 53.4396

DF1 = 1

DF2 = 0

See Also

| |

Statistics and Machine Learning Toolbox Documentation

Support