# coefTest

Class: LinearMixedModel

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

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

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

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

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

`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.4112e-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.4112e-12 ```
```F = 53.0729 ```
```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.0729 ```
```DF1 = 1 ```
```DF2 = 0 ```