lmeVarComp (version 1.1)

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

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.

References

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

Examples

# NOT RUN {
# 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)
# }