rlr.test: Restricted Likelihood Ratio Test and Generalized F-test for...

In lmeVarComp: Testing for a Subset of Variance Components in Linear Mixed Models

Description

rlr.test tests whether certain variance components are zeros using restricted likelihood ratio test and generalized F-test.

Usage

rlr.test(Y, X, Z, Sigma, m0, nsim = 5000L, seed = 130623L)

Arguments

Y	response vector of length n
X	fixed effects design matrix of dimension n by p
Z	a list of random effects design matrices. Each matrix should have n rows.
Sigma	a list of random effects correlation structures. Each matrix should be symmetric and positive definite, and match the dimension of the corresponding random effects design matrix.
m0	an integer indicating the number of nuisance variance components. Should be between 0 and length(Z) - 1. The first m0 variance components will be treated as nuisance.
nsim	number of simulations from the null distribution. If zero, REML estimates are computed but tests are not performed.
seed	a seed to be set before simulating from the null distribution.

Value

A list containing the following components:

RLRT	a vector of the test statistic and the p-value of restricted likelihood ratio test.
GFT	a vector of the test statistic and the p-value of generalized F-test.
H0.estimate	REML estimate of variance components (including the error term) under the null hypothesis.
H1.estimate	REML estimate of variance components (including the error term) under the alternative hypothesis.

Author(s)

Yichi Zhang

References

Zhang, Y., Staicu, A.-M., and Maity, A. (2016). Testing for additivity in non-parametric regression. Canadian Journal of Statistics, 44: 445-462. doi: 10.1002/cjs.11295

Examples

# two-way random effects ANOVA
n1 <- 5L
n2 <- 6L
n0 <- 4L
n <- n1 * n2 * n0
X <- cbind(rep(1, n))
A <- gl(n1, n2 * n0)
Z1 <- model.matrix(~ -1 + A, contrasts.arg = contr.treatment)
B <- rep(gl(n2, n0), n1)
Z2 <- model.matrix(~ -1 + B, contrasts.arg = contr.treatment)
Z3 <- model.matrix(~ -1 + B : A, contrasts.arg = contr.treatment)
set.seed(1L)
Y <- (X %*% 1
  + Z1 %*% rnorm(ncol(Z1), 0, 0.7)
  + Z2 %*% rnorm(ncol(Z2), 0, 0.3)
  + Z3 %*% rnorm(ncol(Z3), 0, 0.5)
  + rnorm(n, 0, 1))
Z <- list(Z1, Z2, Z3)
Sigma <- lapply(Z, function(z) diag(ncol(z)))
# tests interaction effects
rlr.test(Y, X, Z, Sigma, 2L, 2000L, 2L)
# tests overall effects
rlr.test(Y, X, Z, Sigma, 1L, 2000L, 3L)

Example output

$RLRT
stat.obs  p.value 
1.211426 0.116500 

$GFT
stat.obs  p.value 
9.780759 0.118000 

$H0.estimate
[1] 0.35356011 0.05878101 0.85340994

$H1.estimate
[1] 0.34604797 0.04489982 0.08961167 0.79959759

$RLRT
stat.obs  p.value 
3.570463 0.044500 

$GFT
stat.obs  p.value 
17.55381  0.03850 

$H0.estimate
[1] 0.4503803 1.1592004

$H1.estimate
[1] 0.48130011 0.06244882 0.12463621 1.11211868

lmeVarComp documentation built on May 2, 2019, 8:55 a.m.