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

Share:

Description

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

Usage

1
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. (2014). Testing for a Subset of Variance Components in Linear Mixed Models with Application to Testing for Additivity. Submitted.

See Also

pseudo.rlr.test, score.test

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# 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)