Performs permutation tests for random terms in
REML analysis (V.M. Cave).
|Controls printed output (
||Controls the output from the
||What graphs to plot (
||Random term(s) to drop from the full model; no default, must be set|
||Number of permutations to make; default 99|
||Maximum number of extra permutations to make when some
||Seed for random number generation; default 0 continues an existing sequence or, if none, selects a seed automatically|
||Window to use for the graphs; default 3|
||Variates to be analysed|
||Saves the test statistics|
||Saves the p-values|
||Title for the graphs|
||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
VFSTRUCTURE procedures. These define the model settings controlled by the
VSTRUCTURE directives, respectively. The
VAOPTIONS procedure can be used to control some options of the
REML commands used by
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.
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
SAVE parameters. When a single random term is tested,
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,
PROBABILITIES save scalars that store the rLR-based test statistic and its p-value, respectively.
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.
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
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
||prints monitoring information, showing the progress of the analysis.|
||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.|
||prints any error diagnostics from the
VPRINT option controls output from the
REML analyses of the null and alternative models. The settings are the same as those of the
REML directive. By default, nothing is printed.
PLOT option controls graphical output, with settings:
||to plot a histogram of the permuted test statistics, and|
||to produce a kernel density plot of the permuted test statistics.|
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 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 (M1) and reduced model (M0), which omits the random terms specified by
RDROP, are fitted using
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
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
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.
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.
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