Performs permutation tests for random terms in `REML`

analysis (V.M. Cave).

### Options

`PRINT` = string tokens |
Controls printed output (`summary` , `monitoring` , `vdiagnostics` ); default `summ` |
---|---|

`VPRINT` = string tokens |
Controls the output from the `REML` analysis of the full and reduced models (`model` , `components` , `effects` , `means` , `stratumvariances` , `monitoring` , `vcovariance` , `deviance` , `Waldtests` , `missingvalues` , `covariancemodels` ); default `*` i.e. none |

`PLOT` = string tokens |
What graphs to plot (`kerneldensity` , `histogram` ); default `*` |

`MODELDEFINITION` = pointer |
`REML` model definition structure, defined using the `VFMODEL` and `VFSTRUCTURE` procedures, to specify the full model; no default, must be set |

`RDROP` = formula |
Random term(s) to drop from the full model; no default, must be set |

`NTIMES` = scalar |
Number of permutations to make; default 99 |

`NRETRIES` = scalar |
Maximum number of extra permutations to make when some `REML` analyses fail to converge; default `NTIMES` |

`SEED` = scalar |
Seed for random number generation; default 0 continues an existing sequence or, if none, selects a seed automatically |

`WINDOW` = scalar |
Window to use for the graphs; default 3 |

### Parameters

`Y` = variates |
Variates to be analysed |
---|---|

`STATISTICS` = scalars or pointers |
Saves the test statistics |

`PROBABILITIES` = scalars or pointers |
Saves the p-values |

`TITLE` = text |
Title for the graphs |

`SAVE` = pointers |
Saves the test statistics and permuted values |

### Description

`VRPERMTEST`

performs permutation tests to assess whether random terms can be dropped from a linear mixed model, fitted using `REML`

. The procedure implements a pair of permutation tests: one based on the best linear unbiased predictors (BLUP-based), and one based on the residual (or restricted) likelihood ratio test statistic (rLR-based). The rLR-based test enables the simultaneous testing of multiple random terms, whereas the BLUP-based approach can be used only to test the significance of dropping a single random term.

The full model is specified by forming a model-definition structure using the `VFMODEL`

and `VFSTRUCTURE`

procedures. These define the model settings controlled by the `VCOMPONENTS`

and `VSTRUCTURE`

directives, respectively. The `VAOPTIONS`

procedure can be used to control some options of the `REML`

commands used by `VRPERMTEST`

.

The model-definition structure is supplied to `VRPERMTEST`

by the `MODELDEFINITION`

option. The random terms to drop from the full model are defined by a model formula supplied by the `RDROP`

option. If more than one random term is specified, only the rLR-based permutation test is performed.

The `Y`

parameter specifies the variate that is to be modelled. Restrictions and missing values are not allowed in either the y-variate or the explanatory variates or factors. Results can be saved using the `STATISTICS`

, `PROBABILITIES`

and `SAVE`

parameters. When a single random term is tested, `STATISTICS`

and `PROBABILITIES`

save pointers that store the rLR-based and BLUP-based test statistics and their p-values, respectively. Alternatively, when more than one random term is tested, `STATISTICS`

and `PROBABILITIES`

save scalars that store the rLR-based test statistic and its p-value, respectively.

The `SAVE`

parameter can supply a pointer to store the test statistic and its permuted values. The first element of the pointer, indexed by `'permutedT_rLR'`

, is a variate storing the rLR-based test statistic and its permuted values; the first value in the variate is the test statistic from the original data set. When a single random term is tested, a second element, indexed by `'permutedT_BLUP'`

, stores the BLUP-based test statistic (first value) and its permuted values.

The `NTIMES`

option defines how many random permutations to perform; default 99. The `NRETRIES`

option specifies the maximum number of extra samples to take when some `REML`

analyses fail to converge; the default is to use the same number as specified by `NTIMES`

. The `SEED`

option specifies the seed for the random-number generator, used by `RANDOMIZE`

to make the permutations. The default of zero continues the sequence of random numbers from a previous generation or, if this is the first use of the generator in this run of Genstat, it initializes the seed automatically. If you use the same (non-zero) seed more than once, you will get the same random numbers, and hence the same permutations.

Printed output is controlled by the `PRINT`

option, with the following settings.

`monitoring` |
prints monitoring information, showing the progress of the analysis. |
---|---|

`summary` |
prints a summary of the test results: first the seed, the number of permutations and percentage of successful permutations, then a table showing the random term(s) tested, test statistic(s) and the corresponding p-value(s). This is the default. |

`vdiagnostics` |
prints any error diagnostics from the `REML` analyses. |

The `VPRINT`

option controls output from the `REML`

analyses of the null and alternative models. The settings are the same as those of the `PRINT`

option of the `REML`

directive. By default, nothing is printed.

The `PLOT`

option controls graphical output, with settings:

`histogram` |
to plot a histogram of the permuted test statistics, and |
---|---|

`kerneldensity` |
to produce a kernel density plot of the permuted test statistics. |

The `WINDOW`

option defines the window to use for the plots; default 3. By default, nothing is plotted. The `TITLE`

parameter can supply a title for the plots.

Options: `PRINT`

, `VPRINT`

, `PLOT`

, `MODELDEFINITION`

, `RDROP`

, `NTIMES`

, `NRETRIES`

, `SEED`

, `WINDOW`

.

Parameters: `Y`

, `STATISTICS`

, `PROBABILITIES`

, `TITLE`

, `SAVE`

.

### Method

`VRPERMTEST`

is based on the methods of Lee & Braun (2012). Two permutation tests are available. The first test, based on the best linear unbiased predictors (BLUP-based), can be used for inference about a single random effect. The second test, based on the restricted likelihood ratio test statistic (rLR-based), can simultaneously test for the presence of multiple random eﬀects. Both methods involve permuting the weighted marginal residuals. The weights, determined by the Cholesky decomposition of the unit-by-unit variance-covariance matrix, ensure that the marginal residuals are exchangeable under the null hypothesis.

The permutation test proceeds as follows:

1 The full model (*M*_{1}) and reduced model (*M*_{0}), which omits the random terms specified by `RDROP`

, are fitted using `REML`

.

2 The observed test statistic is calculated.

BLUP-based: *T _{BLUP}* = ∑

_{i=}_{1,N}

*b*

_{i}^{2}

where *b*_{1} … *b _{N}* are the estimated BLUPs of the single random term being tested.

rLR-based: *T _{rLR}* = -2 log(

*L*

_{M}_{0}–

*L*

_{M}_{1})

where *L _{M}*

_{0}and

*L*

_{M}_{1}are the restricted likelihoods under the reduced and full models, respectively.

3 The marginal residuals, estimated from the full model, are weighted by (U_{0}′)^{-1}, where U_{0} is the Cholesky decomposition of the unit-by-unit variance-covariance matrix from the reduced model.

4 The weighted errors are permuted using `RANDOMIZE`

.

5 The permuted residuals are unweighted, by multiplication by U_{0}′, and a permuted Y variate (Y*) is obtained.

6 The full and reduced models are refitted with the permuted Y variate, Y*, and permuted values of the test statistics (*T _{BLUP}* and

*T*) are calculated.

_{rLR}7 Steps 4-6 are repeated a maximum of `NTIMES`

+ `NRETRIES`

times.

8 The p-value is given by the proportion of test statistics (including the observed test statistic) greater than the observed test statistic.

The kernel density plot is generated by the `KERNELDENSITY`

procedure, using the method of Sheather & Jones (1991), the default number of grid points, and quantiles calculated at 0.025, 0.25, 0.5, 0.75 and 0.975. The permuted test statistics are plotted using red `+`

symbols along the x-axis, and the location of the test statistic is denoted by a blue line. As the observed test statistic contributes to the null distribution, it is included in the calculation of both the kernel density and histogram.

### Action with `RESTRICT`

Restrictions are not allowed.

### References

Lee, O.E. & Braun, T.M. (2012). Permutation tests for random effects in linear mixed models. *Biometrics*, 68, 486-493.

Sheather, S.J. & Jones, M.C. (1991). A reliable data-based bandwidth selection method for kernel density estimation. *Journal of the Royal Statistical Society, Series B*, 53, 683-690.

### See also

Directives: `REML`

, `VCOMPONENTS`

.

Procedures: `VBOOTSTRAP`

, `VAOPTIONS`

, `VFMODEL`

, `VFSTRUCTURE`

.

Commands for: REML analysis of linear mixed models.

### Example

CAPTION 'VRPERMTEST example',!t('Slate Hall Farm data',\ '(Guide to REML in Genstat, Sections 1.8 and 3.2).');\ STYLE=meta,plain SPLOAD [PRINT=*] '%gendir%/data/slatehall.gsh' " Create model-definition structure with Ar1 (x) Ar1 spatial model, measurement error (plotnumber) and blocking (replicates)." VFMODEL [MODELSTRUCTURE=ModelA; DESCRIPTION='Ar1xAr1 + error + blocking';\ FIXED=variety] fieldrow.fieldcolumn+plotnumber+replicates VFSTRUCTURE [MODELSTRUCTURE=ModelA; TERMS=fieldrow.fieldcolumn]\ 2('AR'); ORDER=1; FACTOR=fieldrow,fieldcolumn " Simultaneous test of plotnumber (measurement error) & replicates (blocking)." VRPERMTEST [PRINT=summary; VPRINT=model,components,deviance;\ PLOT=kerneldensity; MODELDEFINITION=ModelA;\ RDROP=plotnumber+replicates; SEED=16821] yield " Test of the replicates (blocking) random term only " VRPERMTEST [PRINT=summary; VPRINT=model,components,deviance;\ PLOT=kerneldensity; MODELDEFINITION=ModelA;\ RDROP=replicates; SEED=29268] yield " Create model-definition structure without replicates (blocking) random term." VFMODEL [MODELSTRUCTURE=ModelB; IMODELSTRUCTURE=ModelA;\ CHANGEITEMS='random','description';\ DESCRIPTION='Ar1xAr1 + m.error']\ RANDOM=fieldrow.fieldcolumn+plotnumber " Test of plotnumber (measurement error) random term." VRPERMTEST [PRINT=summary; VPRINT=model,components,deviance;\ PLOT=kerneldensity; MODELDEFINITION=ModelB;\ RDROP=plotnumber; SEED=30527] yield