RLRsim (version 2.1-0)

LRTSim: Simulation of the (Restricted) Likelihood Ratio Statistic

Description

These functions simulate values from the (exact) finite sample distribution of the (restricted) likelihood ratio statistic for testing the presence of the variance component (and restrictions of the fixed effects) in a simple linear mixed model with known correlation structure of the random effect and i.i.d. errors. They are usually called by exactLRT or exactRLRT.

Usage

LRTSim(X, Z, q, sqrt.Sigma, seed = NA, nsim = 10000, log.grid.hi = 8,
       log.grid.lo=-10, gridlength=200,
       parallel = c("no", "multicore", "snow"), 
       ncpus = 1L, cl = NULL)
RLRTSim(X, Z, qrX, sqrt.Sigma, lambda0 = NA, seed = NA, nsim = 10000, use.approx=0,
       log.grid.hi=8, log.grid.lo=-10, gridlength=200,
       parallel = c("no", "multicore", "snow"), 
       ncpus = 1L, cl = NULL)

Arguments

X

The fixed effects design matrix of the model under the alternative

Z

The random effects design matrix of the model under the alternative

qrX

The QR-decomposition of the fixed effects design matrix

q

The number of parameters restrictions on the fixed effects (see Details)

sqrt.Sigma

The upper triangular cholesky factor of the correlation matrix of the random effect

lambda0

The value of the ratio of the variance of the random effect and the errors under the null

seed

Specify a seed for set.seed

nsim

Number of values to simulate

use.approx

If 0, the exact distribution is simulated. If between 0 and 1, only the largest eigenvalues whose sum represents at least use.approx*(sum of all eigenvalues) are used.

log.grid.hi

Lower value of the grid on the log scale. See Details

log.grid.lo

Lower value of the grid on the log scale. See Details

gridlength

Length of the grid for the grid search over lambda. See Details

parallel

The type of parallel operation to be used (if any). If missing, the default is "no parallelization").

ncpus

integer: number of processes to be used in parallel operation: typically one would chose this to the number of available CPUs. Defaults to 1, i.e., no parallelization.

cl

An optional parallel or snow cluster for use if parallel = "snow". If not supplied, a cluster on the local machine is created for the duration of the call.

Value

• A vector containig the the simulated values of the (R)LRT under the null, with attribute 'lambda' giving $\arg\min(f(\lambda))$ (see Crainiceanu, Ruppert (2004)) for the simulations.

Details

The model under the alternative must be a linear mixed model $y=X\beta+Zb+\varepsilon$ with a single random effect $b$ with known correlation structure $Sigma$ and i.i.d errors. The simulated distribution of the likelihood ratio statistic was derived by Crainiceanu & Ruppert (2004). The simulation algorithm uses a gridsearch over a log-regular grid of values of $\lambda=\frac{Var(b)}{Var(\varepsilon)}$ to maximize the likelihood under the alternative for nsim realizations of $y$ drawn under the null hypothesis. log.grid.hi and log.grid.lo are the lower and upper limits of this grid on the log scale. gridlength is the number of points on the grid.\ These are just wrapper functions for the underlying C code.

References

Crainiceanu, C. and Ruppert, D. (2004) Likelihood ratio tests in linear mixed models with one variance component, Journal of the Royal Statistical Society: Series B,66,165--185. Scheipl, F. (2007) Testing for nonparametric terms and random effects in structured additive regression. Diploma thesis.\ http://www.statistik.lmu.de/~scheipl/downloads/DIPLOM.zip.

Scheipl, F., Greven, S. and Kuechenhoff, H (2008) Size and power of tests for a zero random effect variance or polynomial regression in additive and linear mixed models, Computational Statistics & Data Analysis, 52(7):3283-3299

See Also

exactLRT, exactRLRT for tests

Examples

library(lme4)
g <- rep(1:10, e = 10)
x <- rnorm(100)
y <- 0.1 * x + rnorm(100)
m <- lmer(y ~ x + (1|g), REML=FALSE)
m0 <- lm(y ~ 1)

(obs.LRT <- 2*(logLik(m)-logLik(m0)))
X <- getME(m,"X")
Z <- t(as.matrix(getME(m,"Zt")))
sim.LRT <- LRTSim(X, Z, 1, diag(10))
(pval <- mean(sim.LRT > obs.LRT))