tests/nbinom.R

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

if (.Platform$OS.type != "windows") {
    library(lme4)
    cat("lme4 testing level: ", testLevel <- lme4:::testLevel(), "\n")


    getNBdisp <- function(x) getME(x,"glmer.nb.theta")
    ## for now, use hidden functions [MM: this is a sign, we should *export* them]
    refitNB   <- lme4:::refitNB

    simfun <- function(sd.u=1, NBtheta=0.5,
                       nblock = 25,
                       fform = ~x,
                       beta = c(1,2),
                       nrep = 40, seed) {
        levelset <- c(LETTERS,letters)
        stopifnot(2 <= nblock, nblock <= length(levelset))
        if (!missing(seed)) set.seed(seed)
        ntot <- nblock*nrep
        d1 <- data.frame(x = runif(ntot),
                         f = factor(rep(levelset[1:nblock], each=nrep)))
        u_f <- rnorm(nblock, sd=sd.u)
        X <- model.matrix(fform, data=d1)
        transform(d1, z = rnbinom(ntot, mu = exp(X %*% beta + u_f[f]), size = NBtheta))
    }

    ##' simplified logLik() so we can compare with "glmmADMB" (and other) results
    logLik.m <- function(x) {
        L <- logLik(x)
        attributes(L) <- attributes(L)[c("class","df","nobs")]
        L
    }

    if (testLevel > 1) withAutoprint({
        set.seed(102)
        d.1 <- simfun()
        t1 <- system.time(g1 <- glmer.nb(z ~ x + (1|f), data=d.1, verbose=TRUE))
        g1
        d1 <- getNBdisp(g1)
        (g1B <- refitNB(g1, theta = d1))
        (ddev <- deviance(g1) - deviance(g1B))
        (reld <- (fixef(g1) - fixef(g1B)) / fixef(g1))
        stopifnot(abs(ddev) < 1e-6,   # was 6.18e-7, 1.045e-6, -6.367e-5, now 0
                  abs(reld) < 1e-6)# 0, then 4.63e-6,  now 0
        ## 2 Aug 2015: ddev==reld==0 on 32-bit Ubuntu 12.04

        if(FALSE) {
            ## comment out to avoid R CMD check warning :
            ## library(glmmADMB)
            t2 <- system.time(g2 <- glmmadmb(z~x+(1|f),
                                             data = d.1, family="nbinom"))
            ## matrix not pos definite in sparse choleski
            t2 # 17.1 sec elapsed
            glmmADMB_vals <- list(fixef= fixef(g2),
                                  LL   = logLik(g2),
                                  theta= g2$alpha)
        } else {
            glmmADMB_vals <-
                list(fixef = c("(Intercept)" = 0.928710, x = 2.05072),
                     LL = structure(-2944.62, class = "logLik", df = 4, nobs = 1000L),
                     theta = 0.4487)
        }


        stopifnot(exprs = {
            all.equal(   d1,         glmmADMB_vals$ theta, tolerance=0.003) #   0.0015907
            all.equal(fixef(g1B),    glmmADMB_vals$ fixef, tolerance=0.02)# was 0.009387 !
            ## Ubuntu 12.04/32-bit: 0.0094
            all.equal(logLik.m(g1B), glmmADMB_vals$ LL,    tolerance=1e-4)# 1.681e-5; Ubuntu 12.04/32-b: 1.61e-5
        })

    })## end if( testLevel > 1 )

    if(FALSE) { ## simulation study --------------------

        ## library(glmmADMB) ## avoid R CMD check warning
        simsumfun <- function(...) {
            d <- simfun(...)
            t1 <- system.time(g1 <- glmer.nb(z~x+(1|f),data=d))
            t2 <- system.time(g2 <- glmmadmb(z~x+(1|f),
                                             data=d,family="nbinom"))
            c(t.glmer=unname(t1["elapsed"]),nevals.glmer=g1$nevals,
              theta.glmer=exp(g1$minimum),
              t.glmmadmb=unname(t2["elapsed"]),theta.glmmadmb=g2$alpha)
        }

        ## library(plyr)
        ## sim50 <- raply(50,simsumfun(),.progress="text")
        save("sim50",file="nbinomsim1.RData")
        ## library(reshape)
        ## m1 <- melt(data.frame(run=seq(nrow(sim50)),sim50),id.var="run")
        ## m1 <- data.frame(m1,colsplit(m1$variable,"\\.",c("v","method")))
        ## m2 <- cast(subset(m1,v=="theta",select=c(run,value,method)),
        ##           run~method)

        library(ggplot2)
        ggplot(subset(m1,v=="theta"),aes(x=method,y=value))+
            geom_boxplot()+geom_point()+geom_hline(yintercept=0.5,colour="red")

        ggplot(subset(m1,v=="theta"),aes(x=method,y=value))+
            stat_summary(fun.data=mean_cl_normal)+
            geom_hline(yintercept=0.5,colour="red")

        ggplot(m2,aes(x=glmer-glmmadmb))+geom_histogram()
        ## glmer is slightly more biased (but maybe the MLE itself is biased???)

    }## end{simulation study}-------------------------

### epilepsy example:
    data(epil, package="MASS")
    epil2 <- transform(epil,
                       Visit  = (period-2.5)/5,
                       Base   = log(base/4),
                       Age    = log(age),
                       subject= factor(subject))

    if(FALSE) {
        ## comment out to avoid R CMD check warning :
        ## library(glmmADMB)
        t3 <- system.time(g3  <- glmmadmb(y~Base*trt+Age+Visit+(Visit|subject),
                                          data=epil2, family="nbinom")) # t3 : 8.67 sec
        glmmADMB_epil_vals <- list(fixef= fixef(g3),
                                   LL   = logLik(g3),
                                   theta= g3$alpha)
    } else {
        glmmADMB_epil_vals <-
            list(fixef =
                     c("(Intercept)"= -1.33, "Base"=0.8839167, "trtprogabide"= -0.9299658,
                       "Age"= 0.4751434, "Visit"=-0.2701603, "Base:trtprogabide"=0.3372421),
                 LL = structure(-624.551, class = "logLik", df = 9, nobs = 236L),
                 theta = 7.4702)
    }

    if (testLevel > 2) withAutoprint({
        ## "too slow" for regular testing -- 49 (MM@lynne: 33, then 26, then 14) seconds:
        (t4 <- system.time(g4 <- glmer.nb(y ~ Base*trt + Age + Visit + (Visit|subject),
                                          data = epil2, verbose=TRUE)))
        ## 1.1-7 : Warning in checkConv().. failed .. with max|grad| = 0.0089 (tol = 0.001, comp. 4)
        ## 1.1-21: 2 Warnings:  max|grad| = 0.00859, then 0.1176 (0.002, comp. 1)

        stopifnot(exprs = {
            all.equal(getNBdisp(g4),   glmmADMB_epil_vals$ theta, tolerance= 0.03) # 0.0019777
            all.equal(fixef    (g4),   glmmADMB_epil_vals$ fixef, tolerance= 0.04) # 0.003731 (0.00374 on U 12.04)
            ## FIXME: even df differ (10 vs 9) !
            ##      all.equal(logLik.m(g4), - glmmADMB_epil_vals$ LL,	tolerance= 0.0) ## was 0.0002
            all.equal(logLik.m(g4), # for now {this is not *the* truth, just our current approximation of it}:
                      structure(-624.48418, class = "logLik", df = 10, nobs = 236L))
        })
    })


    cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons''
} ## skip on windows (for speed)