R Under development (unstable) (2023-07-03 r84633 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > 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) > > proc.time() user system elapsed 0.14 0.06 0.18