Nothing
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)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.