1. Home
  2. VRPERMTEST procedure

VRPERMTEST procedure

Performs permutation tests for random terms in REML analysis (V.M. Cave).


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


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


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.




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 effects. 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 (M1) and reduced model (M0), which omits the random terms specified by RDROP, are fitted using REML.

2    The observed test statistic is calculated.

BLUP-based:                           TBLUP = ∑i=1,N bi2

where b1bN are the estimated BLUPs of the single random term being tested.

rLR-based:                              TrLR = -2 log( LM0LM1 )

where LM0 and LM1 are the restricted likelihoods under the reduced and full models, respectively.

3   The marginal residuals, estimated from the full model, are weighted by (U0′)-1, where U0 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 U0′, 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 (TBLUP and TrLR) are calculated.

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.


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

Commands for: REML analysis of linear mixed models.


CAPTION     'VRPERMTEST example',!t('Slate Hall Farm data',\
            '(Guide to REML in Genstat, Sections 1.8 and 3.2).');\ 
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."
            DESCRIPTION='Ar1xAr1 + m.error']\ 
" Test of plotnumber (measurement error) random term."
VRPERMTEST  [PRINT=summary; VPRINT=model,components,deviance;\
            PLOT=kerneldensity; MODELDEFINITION=ModelB;\ 
            RDROP=plotnumber; SEED=30527] yield
Updated on January 12, 2022

Was this article helpful?