library("testthat")
library("lme4")
context("testing fixed-effect design matrices for full rank")
test_that("lmerRank", {
set.seed(101)
n <- 20
x <- y <- rnorm(n)
d <- data.frame(x,y,
z = rnorm(n),
r = sample(1:5, size=n, replace=TRUE),
y2 = y + c(0.001, rep(0,n-1)))
expect_message(fm <- lmer( z ~ x + y + (1|r), data=d),
"fixed-effect model matrix is .*rank deficient")
## test reconstitution of full parameter vector (with NAs)
expect_equal(names(fixef(fm,add.dropped=TRUE)),
c("(Intercept)","x","y"))
expect_equal(fixef(fm,add.dropped=TRUE)[1:2],
fixef(fm))
expect_equal(nrow(anova(fm)), 1L)
expect_error(lmer( z ~ x + y + (1|r), data=d,
control=lmerControl(check.rankX="stop")),
"rank deficient")
expect_error(lmer( z ~ x + y + (1|r), data=d,
control=lmerControl(check.rankX="ignore")),
"not positive definite")
## should work:
expect_is(lmer( z ~ x + y2 + (1|r), data=d), "lmerMod")
d2 <- expand.grid(a=factor(1:4),b=factor(1:4),rep=1:10)
n <- nrow(d2)
d2 <- transform(d2,r=sample(1:5, size=n, replace=TRUE),
z=rnorm(n))
d2 <- subset(d2,!(a=="4" & b=="4"))
expect_error(lmer( z ~ a*b + (1|r), data=d2,
control=lmerControl(check.rankX="stop")),
"rank deficient")
expect_message(fm <- lmer( z ~ a*b + (1|r), data=d2),
"fixed-effect model matrix is rank deficient")
d2 <- transform(d2, ab=droplevels(interaction(a,b)))
## should work:
expect_is(fm2 <- lmer( z ~ ab + (1|r), data=d2), "lmerMod")
expect_equal(logLik(fm), logLik(fm2))
expect_equal(sum(anova(fm)[, "npar"]), anova(fm2)[, "npar"])
expect_equal(sum(anova(fm)[, "Sum Sq"]), anova(fm2)[, "Sum Sq"])
})
test_that("glmerRank", {
set.seed(111)
n <- 100
x <- y <- rnorm(n)
d <- data.frame(x, y,
z = rbinom(n,size=1,prob=0.5),
r = sample(1:5, size=n, replace=TRUE),
y2 = ## y + c(0.001,rep(0,n-1)), ## too small: get convergence failures
## FIXME: figure out how small a difference will still fail?
rnorm(n))
expect_message(fm <- glmer( z ~ x + y + (1|r), data=d, family=binomial),
"fixed-effect model matrix is rank deficient")
expect_error(glmer( z ~ x + y + (1|r), data=d, family=binomial,
control=glmerControl(check.rankX="stop")),
"rank deficient.*rank.X.")
expect_is(glmer( z ~ x + y2 + (1|r), data=d, family=binomial), "glmerMod")
})
test_that("nlmerRank", {
set.seed(101)
n <- 1000
nblock <- 15
x <- abs(rnorm(n))
y <- rnorm(n)
z <- rnorm(n,mean=x^y)
r <- sample(1:nblock, size=n, replace=TRUE)
d <- data.frame(x,y,z,r)
## save("d","nlmerRank.RData") ## see what's going on with difference in contexts
fModel <- function(a,b) (exp(a)*x)^(b*y)
fModf <- deriv(body(fModel), namevec = c("a","b"),
func = fModel)
fModel2 <- function(a,b,c) (exp(a+c)*x)^(b*y)
fModf2 <- deriv(body(fModel2), namevec = c("a","b","c"),
func = fModel2)
## should be OK: fails in test mode?
nlmer(y ~ fModf(a,b) ~ a|r, d, start = c(a=1,b=1))
## FIXME: this doesn't get caught where I expected
expect_error(nlmer(y ~ fModf2(a,b,c) ~ a|r, d, start = c(a=1,b=1,c=1)),"Downdated VtV")
})
test_that("ranksim", {
set.seed(101)
x <- data.frame(id = factor(sample(10, 100, replace = TRUE)))
x$y <- rnorm(nrow(x))
x$x1 <- 1
x$x2 <- ifelse(x$y<0, rnorm(nrow(x), mean=1), rnorm(nrow(x), mean=-1))
m <- suppressMessages(lmer(y ~ x1 + x2 + (1 | id), data=x))
expect_equal(simulate(m, nsim = 1, use.u = FALSE, newdata=x, seed=101),
simulate(m, nsim = 1, use.u = FALSE, seed=101))
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.