fabian-s/rlrsim

R package for fast & exact (restricted) likelihood ratio tests for mixed and additive models

RLRsim

R package for fast & exact (restricted) likelihood ratio tests for mixed and additive models.

RLRsim implements fast simulation-based exact tests for variance components in mixed and additive models for conditionally Gaussian responses -- i.e., tests for questions like:

• is this random intercept significantly different from 0?
• is this smooth effect significantly nonlinear?
• is this smooth effect significantly different from a constant effect?

The convenience functions exactRLRT and exactLRT can deal with fitted models from packages lme4, nlme, gamm4 and from mgcv's gamm()-function. Workhorse functions LRTSim and RLRTSim accept design matrices as inputs directly and can thus be used more generally to generate exact critical values for the corresponding (restricted) likelihood ratio tests.

The theory behind these tests was first developed in

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.

Power analyses and sensitivity studies for RLRsim can be found in

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.

rlrsim's Issues

test in models with multiple VCs without redoing fit

Needs better formula parsing for nlme and lme4. May mean redoing the interface for exact(R)LRT - deprecate in CRAN version?

installation from source file failed

Dear Fabian,

When I tried to install from source file in my linux machine using
install.packages("PATH/RLRsim_3.0.tar.gz", repos = NULL, type="source")
Error message reported as

*** arch - i386
ERROR: compilation failed for package 'RLRsim'

However, intallation using the windows binary works for my windows machine. What might be the reason?

Thanks a lot.
Ni

change in anova() names breaks

This change breaks exactRLRT: I will try to fork / fix it/ submit a PR. Don't know when a new release of lme4 is happening, but it might be soon-ish.

extract.lmeDesign(lme4_object) not working

Hello,
is it possible to allow the extraction of a s4 class lme4 object?
Currently it gives an error:

library(lme4)
m<-lme4::lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
extract.lmeDesign(m)
Error: $ operator not defined for this S4 class

interface for mgcv fits

just do it already!

model.matrix warning about contrasts

### R code from vignette source 'Ch-LME.Rnw'

###################################################
### code chunk number 1: setup
###################################################
library("MVA")
set.seed(280875)


###################################################
### code chunk number 2: setup
###################################################
library("nlme")

###################################################
### code chunk number 6: ch:LME:timber:tab
###################################################
"timber" <-
matrix(c(0., 0., 0., 0., 0., 0., 0., 0., 2.3799999999999999, 2.6899999999999999, 2.8500000000000001, 2.46,
        2.9700000000000002, 3.96, 3.1699999999999999, 3.3599999999999999, 4.3399999999999999, 4.75,
        4.8899999999999997, 4.2800000000000002, 4.6799999999999997, 6.46, 5.3300000000000001, 5.4500000000000002,
        6.6399999999999997, 7.04, 6.6100000000000003, 5.8799999999999999, 6.6600000000000001, 8.1400000000000006,
        7.1399999999999997, 7.0800000000000001, 8.0500000000000007, 9.1999999999999993, 8.0899999999999999,
        7.4299999999999997, 8.1099999999999994, 9.3499999999999996, 8.2899999999999991, 8.3200000000000003,
        9.7799999999999994, 10.94, 9.7200000000000006, 8.3200000000000003, 9.6400000000000006, 10.720000000000001,
        9.8599999999999994, 9.9100000000000001, 10.970000000000001, 12.23, 11.029999999999999, 9.9199999999999999,
        11.06, 11.84, 11.07, 11.06, 12.050000000000001, 13.19, 12.140000000000001, 11.1, 12.25, 12.85,
        12.130000000000001, 12.210000000000001, 12.98, 14.08, 13.18, 12.23, 13.35, 13.83, 13.15, 13.16, 13.94,
        14.66, 14.119999999999999, 13.24, 14.539999999999999, 14.85, 14.09, 14.050000000000001, 14.74,
        15.369999999999999, 15.09, 14.19, 15.529999999999999, 15.789999999999999, 15.109999999999999,
        14.960000000000001, 16.129999999999999, 16.890000000000001, 16.68, 16.07, 17.379999999999999,
        17.390000000000001, 16.690000000000001, 16.239999999999998, 17.98, 17.780000000000001, 17.940000000000001,
        17.43, 18.760000000000002, 18.440000000000001, 17.690000000000001, 17.34, 19.52, 18.41, 18.219999999999999,
        18.359999999999999, 19.809999999999999, 19.460000000000001, 18.710000000000001, 18.23, 19.969999999999999,
        18.969999999999999, 19.399999999999999, 18.93, 20.620000000000001, 20.050000000000001, 19.539999999999999,
        18.870000000000001)
, nrow = 8, ncol = 15)

slippage <- c((0:10)/10, seq(from = 1.2, to = 1.8, by = 0.2))

colnames(timber) <- paste("s", slippage, sep = "")
timber <- as.data.frame(timber)
timber$specimen <- factor(paste("spec", 1:nrow(timber), sep = ""))

timber.dat <- reshape(timber, direction = "long", idvar = "specimen",
        varying = matrix(colnames(timber)[1:15], nr = 1),
        timevar = "slippage")
names(timber.dat)[3] <- "loads"
timber.dat$slippage <- slippage[timber.dat$slippage]
timber <- timber.dat

###################################################
### code chunk number 11: ch:LME:timber:lme
###################################################
timber.lme <- lme(loads ~ slippage, 
                  random = ~1 | specimen,
                  data = timber, method = "ML")

###################################################
### code chunk number 12: ch:LME:timber:LRT
###################################################
library("RLRsim")

exactRLRT(timber.lme)
#Using restricted likelihood evaluated at ML estimators.
#Refit with method="REML" for exact results.

#	simulated finite sample distribution of RLRT.
	
#	(p-value based on 10000 simulated values)

#data:  
#RLRT = 3.2183e-07, p-value = 0.4419

#Warning message:
#In model.matrix.default(~m$groups[[n.levels - i + 1]] - 1, contrasts.arg = c("contr.treatment",  :
#  non-list contrasts argument ignored

sample of exactRLRT

exactRLRT(m0)$sample is a function instead of simulated test statistics. How should I fix it?