coefTest

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

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

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

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 H₀: Hβ = 0, where β is the fixed-effects vector.

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 H₀: Hβ = C, where β is the fixed-effects vector.

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 H₀: Hβ + KB = C, where β is the fixed-effects vector and B is the random-effects vector.

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

Examples

Test Fixed-Effects Coefficients for Categorical Data

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       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 Tests for Fixed-Effects Coefficients

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

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.

Hypothesis Tests for Fixed- and Random-Effects Coefficients

Load the sample data and convert it to table format.

load flu
flu = dataset2table(flu)

flu=52×11 table
         Date          NE      MidAtl    ENCentral    WNCentral    SAtl     ESCentral    WSCentral     Mtn      Pac     WtdILI
    ______________    _____    ______    _________    _________    _____    _________    _________    _____    _____    ______

    {'10/9/2005' }     0.97    1.025       1.232        1.286      1.082      1.457          1.1      0.981    0.971    1.182 
    {'10/16/2005'}    1.136     1.06       1.228        1.286      1.146      1.644        1.123      0.976    0.917     1.22 
    {'10/23/2005'}    1.135    1.172       1.278        1.536      1.274      1.556        1.236      1.102    0.895     1.31 
    {'10/30/2005'}     1.52    1.489       1.576        1.794       1.59      2.252        1.612      1.321    1.082    1.343 
    {'11/6/2005' }    1.365    1.394        1.53        1.825       1.62      2.059        1.471      1.453    1.118    1.586 
    {'11/13/2005'}     1.39    1.477       1.506          1.9      1.683      1.813        1.464      1.388    1.204     1.47 
    {'11/20/2005'}    1.212    1.231       1.295        1.495      1.347      1.794        1.303      1.371    1.137    1.611 
    {'11/27/2005'}    1.477    1.546       1.557        1.855      1.678      2.159        1.739      1.628    1.443    1.827 
    {'12/4/2005' }    1.285     1.43       1.482        1.635      1.577      1.903         1.53      1.701    1.516    1.776 
    {'12/11/2005'}    1.354     1.45        1.46        1.794      1.583      1.894        1.831      2.364    2.094    1.941 
    {'12/18/2005'}    1.502    1.622       1.638        1.988      1.947       2.22        2.577       3.89     2.66     2.34 
    {'12/25/2005'}     1.86    1.915       1.955         2.38      2.343      3.027        3.219      4.862    2.595    3.086 
    {'1/1/2006'  }    2.114    2.174       2.065        2.557      2.275      2.498        2.644      3.352    2.181     3.26 
    {'1/8/2006'  }    1.815    1.932       1.822        2.046      1.969      1.805        2.189      2.132    1.717    2.613 
    {'1/15/2006' }    1.541    1.695       1.581        2.008      1.718      1.662        2.156      1.694    1.351    2.247 
    {'1/22/2006' }    1.632    1.758       1.711        2.217      1.866      2.194        2.268      1.826    1.384    2.352 
      ⋮

flu contains 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 table. 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 table, 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,NewDataVariableName="FluRate",...
    IndexVariableName="Region");
flu2.Date = categorical(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     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')

pVal = 
1.4770e-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.4770e-12

F = 
52.9730

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

DF1 = 
1

DF2 = 
0

Input Arguments

lme — Linear mixed-effects model

LinearMixedModel object

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

H — Fixed-effects contrasts

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

Example: pVal = coefTest(lme,H)

Data Types: single | double

C — Hypothesized value

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

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: pVal = coefTest(lme,H,C,'DFMethod','satterthwaite')

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

REContrast — Random-effects contrasts

m-by-q matrix

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

pVal — p-value

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

F — F-statistic

scalar value

F-statistic, returned as a scalar value.

DF1 — Numerator degrees of freedom for F

scalar value

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

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

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

DF2 — Denominator degrees of freedom for F

scalar value

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

Version History

Introduced in R2013b