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 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 b1 … bN are the estimated BLUPs of the single random term being tested.
rLR-based: TrLR = -2 log( LM0 – LM1 )
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.
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, VFLC, VFMODEL, VFSTRUCTURE, VPERMTEST.
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