tests/testthat/test-glmer.R

In lme4/lme4: Linear Mixed-Effects Models using 'Eigen' and S4

library("testthat")
library("lme4")
source(system.file("testdata", "lme-tst-funs.R", package="lme4", mustWork=TRUE))# -> uc()
testLevel <- lme4:::testLevel()

gives_error_or_warning <- function (regexp = NULL, all = FALSE, ...)
{
    function(expr) {
        res <- try(evaluate_promise(expr),silent=TRUE)
        no_error <- !inherits(res, "try-error")
        if (no_error) {
            warnings <- res$warnings

            if (!is.null(regexp) && length(warnings) > 0) {
                return(matches(regexp, all = FALSE, ...)(warnings))
            } else {
                return(expectation(length(warnings) > 0, "no warnings or errors given",
                            paste0(length(warnings), " warnings created")))
            }
        }
        if (!is.null(regexp)) {
            return(matches(regexp, ...)(res))
        }
        else {
            expectation(TRUE, "no error thrown", "threw an error")
        }
    }
}
    ## expect_that(stop("foo"),gives_error_or_warning("foo"))
    ## expect_that(warning("foo"),gives_error_or_warning("foo"))
    ## expect_that(TRUE,gives_error_or_warning("foo"))
    ## expect_that(stop("bar"),gives_error_or_warning("foo"))
    ## expect_that(warning("bar"),gives_error_or_warning("foo"))

if(testLevel > 1) {
context("fitting glmer models")
test_that("glmer", {
    set.seed(101)
    d <- data.frame(z=rbinom(200,size=1,prob=0.5),
                    f=factor(sample(1:10,200,replace=TRUE)))
    ## Using 'method=*' defunct in 2019-05 (after 6 years of deprecation)
    ## expect_warning(glmer(z~ 1|f, d, family=binomial, method="abc"),"Use the nAGQ argument")
    ## expect_warning(glmer(z~ 1|f, d, family=binomial, method="Laplace"),"Use the nAGQ argument")
    ##sp expect_warning(glmer(z~ 1|f, d, sparseX=TRUE),"has no effect at present")
    expect_that(gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                             data = cbpp, family = binomial), is_a("glmerMod"))
    expect_that(gm1@resp,                               is_a("glmResp"))
    expect_that(gm1@pp,                                 is_a("merPredD"))
    expect_equal(ge1 <- unname(fixef(gm1)),             c(-1.39854982537216, -0.992335519118859,
                                                          -1.12867532780426, -1.58030423764517),
                                                         tolerance=5e-4)
    expect_equal(c(VarCorr(gm1)[[1]]),                  0.41245527438386, tolerance=6e-4)
### expect_that(family(gm1),                            equals(binomial()))
### ?? binomial() has an 'initialize' component ... and the order is different
    expect_equal(deviance(gm1),                         73.47428, tolerance=1e-5)
    ## was -2L = 184.05267459802
    expect_equal(sigma(gm1),                            1)
    expect_equal(extractAIC(gm1),                       c(5, 194.052674598026), tolerance=1e-5)

    expect_equal(theta <- unname(getME(gm1, "theta")),  0.642226809144453, tolerance=6e-4)
    expect_that(X  <- getME(gm1, "X"),                  is_equivalent_to(
        model.matrix(model.frame(~ period, data=cbpp), cbpp)))
    expect_that(Zt <- getME(gm1, "Zt"),                 is_a("dgCMatrix"))
    expect_equal(dim(Zt),                               c(15L, 56L))
    expect_equal(Zt@x,                                  rep.int(1, 56L))
    expect_that(Lambdat <- getME(gm1, "Lambdat"),       is_a("dgCMatrix"))
    expect_equivalent(as(Lambdat, "matrix"),            diag(theta, 15L, 15L))

    expect_is(gm1_probit <- update(gm1,family=binomial(link="probit")),"merMod")
    expect_equal(family(gm1_probit)$link,"probit")

    ## FIXME: test user-specified/custom family?

    expect_error(glFormula(cbind(incidence, size - incidence) ~ period + (1 | herd),
                             data = subset(cbpp, herd==levels(herd)[1]), family = binomial),
                 "must have > 1")
    expect_warning(glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                         data = subset(cbpp, herd %in% levels(herd)[1:4]),
                         family = binomial,
                         control=glmerControl(check.nlev.gtreq.5="warning")),
                   "< 5 sampled levels")
    expect_warning(fm1. <- glmer(Reaction ~ Days + (Days|Subject), sleepstudy),
                   regexp="calling .* with family=gaussian .* as a shortcut")
    options(warn=2)
    options(glmerControl=list(junk=1,check.conv.grad="ignore"))
    expect_warning(glmer(z~ 1|f, d, family=binomial),
                   "some options")
    options(glmerControl=NULL)
    cbppX <- transform(cbpp,prop=incidence/size)
    expect_is(glmer(prop ~ period + (1 | herd),
		      data = cbppX, family = binomial, weights=size),
	      "glmerMod")
    expect_is(glmer(prop ~ period + (1 | herd),
		      data = cbppX, family = binomial, weights=size, start=NULL),
	      "glmerMod")
    expect_is(glmer(prop ~ period + (1 | herd),
		      data = cbppX, family = binomial, weights=size, verbose=0L),
	      "glmerMod")
    expect_is(glmer(prop ~ period + (1 | herd),
		      data = cbppX, family = binomial, weights=size, subset=TRUE),
	      "glmerMod")
    expect_is(glmer(prop ~ period + (1 | herd),
		      data = cbppX, family = binomial, weights=size, na.action="na.exclude"),
	      "glmerMod")

    expect_is(glmer(prop ~ period + (1 | herd),
    		      data = cbppX, family = binomial, weights=size, offset=rep(0,nrow(cbppX))),
    "glmerMod")
    expect_is(glmer(prop ~ period + (1 | herd),
		      data = cbppX, family = binomial, weights=size, contrasts=NULL),
	      "glmerMod")
    expect_is(glmer(prop ~ period + (1 | herd),
		      data = cbppX, family = binomial, weights=size, devFunOnly=FALSE),
	      "glmerMod")
    expect_is(glmer(prop ~ period + (1 | herd),
		      data = cbppX, family = binomial, weights=size,
		    control=glmerControl(optimizer="Nelder_Mead")),
	      "glmerMod")
    expect_is(glmer(prop ~ period + (1 | herd),
		      data = cbppX, family = binomial, weights=size, control=glmerControl()),
	      "glmerMod")
    options(warn=0)
    expect_error(glmer(prop ~ period + (1 | herd),
                       data = cbppX, family = binomial, weights=size, junkArg=TRUE),
                 "unused argument")
if(FALSE) { ## Hadley broke this
    expect_warning(glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                         data = cbpp, family = binomial,
                         control=list()),
                   "instead of passing a list of class")
    expect_warning(glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                         data = cbpp, family = binomial,
                         control=lmerControl()),
                   "instead of passing a list of class")
}
    ##
    load(system.file("testdata","radinger_dat.RData",package="lme4"))
    mod <- glmer(presabs~predictor+(1|species),family=binomial,
                 radinger_dat)
    expect_is(mod,"merMod")
    ## tolerance: 32-bit Windows (CRAN) reported ave.diff of 5.33e-8
    ## 64-bit Win-builder r73242 now reports ave. diff of 1.31e-5 ...
    expect_equal(unname(fixef(mod)), c(0.5425528,6.4289962), tolerance = 1e-4)
    set.seed(101)
    ## complete separation case
    d <- data.frame(y=rbinom(1000,size=1,p=0.5),
                    x=runif(1000),
                    f=factor(rep(1:20,each=50)),
                    x2=rep(0:1,c(999,1)))
    expect_message(mod2 <- glmer(y~x+x2+(1|f),data=d,family=binomial,
                                 control=glmerControl(check.conv.hess="ignore",
                                                      check.conv.grad="ignore")),
                   "singular")
    expect_equal(unname(fixef(mod2))[1:2],
                 c(-0.10036244,0.03548523), tolerance=1e-4)
    expect_true(unname(fixef(mod2)[3] < -10))
    expect_message(mod3 <- update(mod2, family=binomial(link="probit")),
                   "singular")
    # singular Hessian warning
    expect_equal(unname(fixef(mod3))[1:2], c(-0.062889, 0.022241), tolerance=1e-4)
    expect_true(fixef(mod3)[3] < -4)
    expect_message(mod4 <- update(mod2, family=binomial(link="cauchit"),
                   control=glmerControl(check.conv.hess="ignore",
                                        check.conv.grad="ignore")))#--> singular Hessian warning

    ## on-the-fly creation of index variables
    if (FALSE) {
        ## FIXME: fails in testthat context -- 'd' is not found
        ##  in the parent environment of glmer() -- but works fine
        ## otherwise ...
        set.seed(101)
        d <- data.frame(y1=rpois(100,1),  x=rnorm(100), ID=1:100)
        fit1 <- glmer(y1 ~ x+(1|ID),data=d,family=poisson)
        fit2 <- update(fit1, .~ x+(1|rownames(d)))
        expect_equal(unname(unlist(VarCorr(fit1))),
                     unname(unlist(VarCorr(fit2))))
    }

    ##
    if(testLevel > 2) {
        load(system.file("testdata","mastitis.rda",package="lme4"))
        t1 <- system.time(g1 <-
                              suppressWarnings(glmer(NCM ~ birth + calvingYear + (1|sire) + (1|herd),
                                mastitis, poisson,
                                ## current (2014-04-24) default: --> Warning
                                control=glmerControl(  # max|grad| = 0.021 ..
                                    optimizer=c("bobyqa","Nelder_Mead")))))

        t2 <- system.time(g2 <- suppressWarnings(update(g1,
                         control=glmerControl(optimizer="bobyqa"))))
        ## rbind(t1,t2)[,"elapsed"]
        ## 20 (then 13.0) seconds N-M vs 8 (then 4.8) seconds bobyqa ...
        ## print(t1[3] / t2[3]) # 0.37; => 1.25 should be on the safe side
        expect_lte(t2[3], 1.25 * t1[3])
        ## problem is fairly ill-conditioned so parameters
        ##  are relatively far apart even though likelihoods are OK
        expect_equal(logLik(g1),logLik(g2),tolerance=2e-7)
    }
    ## test bootstrap/refit with nAGQ>1
    gm1AGQ <- update(gm1,nAGQ=2)
    s1 <- simulate(gm1AGQ)
    expect_equal(attr(bootMer(gm1AGQ,fixef),"bootFail"),0)

    ## do.call(new,...) bug
    new <- "foo"
    expect_that(gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                             data = cbpp, family = binomial), is_a("glmerMod"))
    rm("new")

    ## test issue #47, from Wolfgang Viechtbauer
    ## create some data
    n <- 100
    ai <- rep(0:1, each = n/2)
    bi <- 1-ai
    ci <- c(rep(0,42), rep(1,8), rep(0,18), rep(1,32))
    di <- 1-ci
    event <- c(rbind(ai,ci))
    group <- rep(c(1,0), times=n)
    id    <- rep(1:n, each=2)
    gm3 <- glmer(event ~ group + (1 | id), family=binomial, nAGQ=21)
    sd3 <- sqrt(diag(vcov(gm3)))
    expect_equal(uc(`(Intercept)` = 0.42542855, group = 0.42492581), sd3, tolerance=1e-5) # 7e-9 {Lnx}
    expect_warning(vcov(gm3,use.hessian=FALSE), "finite-difference Hessian")
    expect_equal(suppressWarnings(sqrt(diag(vcov(gm3,use.hessian=FALSE)))),
                 uc(`(Intercept)` = 0.3840921, group = 0.3768747), tolerance=1e-7) # 6.5e-8
    expect_equal(sd3, unn(coef(summary(gm3))[,"Std. Error"]))
    ## test non-pos-def finite-difference Hessian ...
    if(getRversion() > "3.0.0") {
        ## saved fits are not safe with old R versions
        L <- load(system.file("testdata","polytomous_vcov_ex.RData",
                              package="lme4", mustWork=TRUE))
        expect_warning(vcov(polytomous_vcov_ex),"falling back to var-cov")
    }

    ## damage Hessian to make it singular
    ## (example thanks to J. Dushoff)
    gm1H <- gm1
    gm1H@optinfo$derivs$Hessian[5,] <- 0
    expect_warning(vcov(gm1H),"falling back to var-cov")

    ## test convergence warnings
    L <- load(system.file("testdata","gopherdat2.RData",
                          package="lme4", mustWork=TRUE))
    g0 <- glmer(shells~prev + (1|Site)+offset(log(Area)),
                family=poisson, data=Gdat)
    ## fit year as factor: OK
    gc <- glmerControl(check.conv.grad="stop")
    expect_is(update(g0,.~.+factor(year), control=gc), "glmerMod")
    ## error/warning with year as numeric:
    ## don't have full knowledge of which platforms lead to which
    ##     results, and can't detect whether we're running on valgrind,
    ##     which changes the result on 32-bit linux ...
    ## SEGFAULT on MacOS? why?
    if (FALSE) {
       expect_that(update(g0,.~.+year),
         gives_error_or_warning("(failed to converge|pwrssUpdate did not converge)"))
    }
    ## ("(failed to converge|pwrssUpdate did not converge in)"))
    ## if (sessionInfo()$platform=="i686-pc-linux-gnu (32-bit)") {
    ##     expect_warning(update(g0, .~. +year), "failed to converge")
    ## } else {
    ##     ## MacOS x86_64-apple-darwin10.8.0 (64-bit)
    ##     ## MM's platform
    ##     ## "pwrssUpdate did not converge in (maxit) iterations"
    ##     expect_error(update(g0, .~. +year), "pwrssUpdate did not converge in")
    ## }
    ## OK if we scale & center it
    expect_is(update(g0,.~. + scale(year), control=gc), "glmerMod")
    ## not OK if we scale and don't center
    expect_warning(update(g0,.~. + scale(year,center=FALSE)),
                   "failed to converge with max|grad|")
    ## OK if center and don't scale
    expect_is(update(g0,.~. + scale(year,center=TRUE,scale=FALSE),
                     control=gc),
              "glmerMod")

    ## try higher-order AGQ
    expect_is   (update(gm1,nAGQ=90), "glmerMod")
    expect_error(update(gm1,nAGQ=101),"ord < 101L")

    ## non-numeric response variables
    ss <- transform(sleepstudy, Reaction = as.character(Reaction))
    expect_error(glmer(Reaction~(1|Days), family="poisson", data=ss),
                 "response must be numeric")
    expect_error(glmer(Reaction~(1|Days), family="binomial", data=ss),
                 "response must be numeric or factor")
    ss2 <- transform(ss,rr=rep(c(TRUE,FALSE),length.out=nrow(ss)))
    ## should work OK with logical too
    expect_is(glmer(rr~(1|Days),family="binomial",data=ss2),"merMod")

    ## starting values with log(.) link -- thanks to Eric Weese @ Yale:
    grp <- rep(letters[1:5], 20); set.seed(1); x <- rnorm(100)
    expect_error(glmer(x ~ 1 + (1|grp), family=gaussian(link="log")),
		 "valid starting values")

    ## related to GH 231
    ## fails on some platforms, skip for now
    if (FALSE) {
        rr <- gm1@resp$copy()
        ff <- setdiff(ls(gm1@resp),c("copy","initialize","initialize#lmResp","ptr",
                                     "updateMu","updateWts","resDev","setOffset","wrss"))
        for (i in ff) {
            expect_equal(gm1@resp[[i]],rr[[i]])
        }
    }

    ## bad start case
    load(system.file("testdata","fakesim.RData",package="lme4"))
    rfit <- glmer(Inew/S ~ R0-1 + offset(log(I/N)) + (1|R0:trial)
        , family=binomial(link="cloglog")
	, data=dat
	, weight=S
	, control=glmerControl(optimizer="bobyqa",
                               nAGQ0initStep=FALSE)
	, start = list(fixef=c(0,0,0),theta=1))
    expect_equal(exp(fixef(rfit)),
                 c(R01 = 1.2735051, R02 = 2.0330635, R03 = 2.9764088),
                 tolerance=1e-5)

    ## gaussian with log link and zero values in response ...
    ##  fixed simulation code, passing mustart properly
    dd <- expand.grid(x = seq(-2,3,length.out=10),
                      f = factor(1:10))
    dd$y <- simulate(~x+(1|f),
                 family=gaussian(link="log"),
                 newdata=dd,
                 newparams=list(beta=c(0,1),theta=1,sigma=1),
                 seed=101)[[1]]
    dd$y <- pmax(dd$y,0)
    expect_error(glmer (y ~ x + (1|f),
                        family = gaussian(link="log"),
                        data=dd),"cannot find valid starting values")
    g1 <- glmer (y ~ x + (1|f),
                 family = gaussian(link="log"),
                 data=dd,
                 mustart=pmax(dd$y,0.1))
    msum <- c(fixef(g1),unlist(c(VarCorr(g1))),c(logLik(g1)))
    expect_equal(msum,
                 c(`(Intercept)` = 0.23389405, x = 1.0017436, f = 0.24602992, -156.7773),
                 tolerance=1e-5)

    ## GH 415
    expect_warning(glmer (round(Reaction) ~ Days + (1|Subject),
                          data=sleepstudy[1:100,],
                          family=poisson,
                          control=lmerControl()),
                   "please use glmerControl")

})
} ## testlevel>1

test_that("glmer with etastart",
{
    ## make sure etastart is passed through
    ## (fixed in commit b6fb1ac83885ff06 but never tested?)
    m1 <- glmer(incidence/size ~ period + (1|herd),
                weights = size,
                family = binomial,
                data = cbpp)
    m1E <- update(m1, etastart = rep(1, nrow(cbpp)))
    expect_true(!identical(fixef(m1), fixef(m1E)))
})