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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.