[Lme4-commits] r1630 - in pkg/lme4Eigen: R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 29 12:57:05 CET 2012
Author: mmaechler
Date: 2012-02-29 12:57:05 +0100 (Wed, 29 Feb 2012)
New Revision: 1630
Modified:
pkg/lme4Eigen/R/bootMer.R
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/R/nbinom.R
pkg/lme4Eigen/R/profile.R
pkg/lme4Eigen/man/mcmcsamp.Rd
pkg/lme4Eigen/tests/getME.R
pkg/lme4Eigen/tests/nbinom.R
Log:
a typo in tests/getME.R .. comments and more checking in tests/nbinom.R
Modified: pkg/lme4Eigen/R/bootMer.R
===================================================================
--- pkg/lme4Eigen/R/bootMer.R 2012-02-28 23:26:26 UTC (rev 1629)
+++ pkg/lme4Eigen/R/bootMer.R 2012-02-29 11:57:05 UTC (rev 1630)
@@ -110,7 +110,7 @@
stop("bootMer currently only handles functions that return numeric vectors")
mle <- list(beta = getME(x,"beta"), theta = getME(x,"theta"))
- if (lme4Eigen:::isLMM(x)) mle <- c(mle,list(sigma = sigma(x)))
+ if (isLMM(x)) mle <- c(mle,list(sigma = sigma(x)))
## FIXME: what about GLMMs with scale parameters??
## FIXME: remove prefix when incorporated in package
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2012-02-28 23:26:26 UTC (rev 1629)
+++ pkg/lme4Eigen/R/lmer.R 2012-02-29 11:57:05 UTC (rev 1630)
@@ -1187,21 +1187,23 @@
}
##' @S3method refit merMod
-refit.merMod <- function(object, newresp, ...) {
- rr <- object at resp$copy()
+refit.merMod <- function(object, newresp= model.response(model.frame(object)),
+ ...)
+{
+ rr <- object at resp$copy()
if (isGLMM(object) && rr$family$family=="binomial") {
- if (is.matrix(newresp) && ncol(newresp)==2) {
- ntot <- rowSums(newresp)
- ## FIXME: test what happens for (0,0) columns
- newresp <- newresp[,1]/ntot
- rr$setWeights(ntot)
- }
- if (is.factor(newresp)) {
- ## FIXME: would be better to do this consistently with
- ## whatever machinery is used in glm/glm.fit/glmer ... ??
- newresp <- as.numeric(newresp)-1
- }
+ if (is.matrix(newresp) && ncol(newresp)==2) {
+ ntot <- rowSums(newresp)
+ ## FIXME: test what happens for (0,0) columns
+ newresp <- newresp[,1]/ntot
+ rr$setWeights(ntot)
+ }
+ if (is.factor(newresp)) {
+ ## FIXME: would be better to do this consistently with
+ ## whatever machinery is used in glm/glm.fit/glmer ... ??
+ newresp <- as.numeric(newresp)-1
+ }
}
stopifnot(length(newresp <- as.numeric(as.vector(newresp))) == length(rr$y))
rr$setResp(newresp)
@@ -1212,8 +1214,8 @@
devlist <- list(pp=pp, resp=rr, u0=pp$u0, beta0=pp$beta0, verbose=0L,
dpars=seq_len(nth))
if (isGLMM(object)) {
- devlist <- c(devlist,tolPwrss=unname(dc$cmp["tolPwrss"]),
- nAGQ=unname(nAGQ))
+ devlist <- c(devlist,tolPwrss=unname(dc$cmp["tolPwrss"]),
+ nAGQ=unname(nAGQ))
}
ff <- mkdevfun(list2env(devlist),nAGQ=nAGQ)
xst <- rep.int(0.1, nth)
Modified: pkg/lme4Eigen/R/nbinom.R
===================================================================
--- pkg/lme4Eigen/R/nbinom.R 2012-02-28 23:26:26 UTC (rev 1629)
+++ pkg/lme4Eigen/R/nbinom.R 2012-02-29 11:57:05 UTC (rev 1630)
@@ -2,6 +2,7 @@
##' @importFrom MASS theta.ml
## should be getME(object,"NBdisp") ?
+## MM: should the *user* use it? if yes, consider method sigma() or AIC() ?
getNBdisp <- function(object) {
get(".Theta",envir=environment(object at resp$family$aic))
}
@@ -13,39 +14,33 @@
ff <- setdiff(names(getRefClass("glmResp")$fields()),c("Ptr","family"))
arg1 <- lapply(ff,object at resp$field)
names(arg1) <- ff
- newresp <- do.call(glmResp$new,c(arg1,
- list(family=negative.binomial(theta=theta))))
+ newresp <- do.call(glmResp$new,
+ c(arg1, list(family=negative.binomial(theta=theta))))
object at resp <- newresp
object
}
refitNB <- function(object,theta) {
- orig_theta <- getNBdisp(object)
- object <- setNBdisp(object,theta) ## new/copied object
- refit(object,newresp=model.response(model.frame(object)))
- ## FIXME: should refit() take this response as a default??
- ## Yes, I think that is a good idea. DB 2012-02-23
+ refit(setNBdisp(object,theta))
}
optTheta <- function(object,
interval=c(-5,5),
maxit=20,
- debug=FALSE) {
+ verbose=FALSE) {
lastfit <- object
evalcnt <- 0
optval <- optimize(function(t) {
- ## FIXME: kluge to retain last value and evaluation count
+ ## FIXME: kluge to retain last value and evaluation count
## Perhaps use a reference class object to keep track of this
## auxilliary information? DB
- L <- -logLik(lastfit <<- refitNB(lastfit,theta=exp(t)))
- evalcnt <<- evalcnt+1
- if (debug) {
- cat(evalcnt,exp(t),L,"\n")
- }
+ L <- -logLik(lastfit <<- refitNB(lastfit,theta=exp(t)))
+ evalcnt <<- evalcnt+1
+ if (verbose) cat(evalcnt,exp(t),L,"\n")
L
- },interval=interval)
+ }, interval=interval)
stopifnot(all.equal(optval$minimum,log(getNBdisp(lastfit))))
- ## FIXME: return eval count info somewhere else?
+ ## FIXME: return eval count info somewhere else? MM: new slot there, why not?
attr(lastfit,"nevals") <- evalcnt
lastfit
}
@@ -60,16 +55,22 @@
trace = control$trace > 2)
}
-## wrapper for glmer stuff
+## FIXME: really should document glmer.nb() on the same help page as glmer()
+## I (MM) don't know how to use roxygen for that..
+
+##' glmer() for Negative Binomial
+##' @param ... formula, data, etc: the arguments for \code{\link{glmer}(..)} (apart from \code{family}!).
+##' @param interval interval in which to start the optimization
+##' @param verbose logical indicating how much progress information should be printed.
##' @export
-glmer.nb <- function(...,
- interval=NULL,
- debug=FALSE) {
- g0 <- glmer(...,family=poisson)
+glmer.nb <- function(..., interval = log(th)+c(-3,3), verbose=FALSE)
+{
+ g0 <- glmer(..., family=poisson)
th <- est_theta(g0)
- g1 <- update(g0,family=negative.binomial(theta=th))
- if (is.null(interval)) interval <- log(th)+c(-3,3)
- optTheta(g1,interval=interval,debug=debug)
+ if(verbose) cat("th := est_theta(glmer(..)) =", format(th),"\n")
+ g1 <- update(g0, family = negative.binomial(theta=th))
+ ## if (is.null(interval)) interval <- log(th)+c(-3,3)
+ optTheta(g1, interval=interval, verbose=verbose)
}
## do we want to facilitate profiling on theta??
Modified: pkg/lme4Eigen/R/profile.R
===================================================================
--- pkg/lme4Eigen/R/profile.R 2012-02-28 23:26:26 UTC (rev 1629)
+++ pkg/lme4Eigen/R/profile.R 2012-02-29 11:57:05 UTC (rev 1630)
@@ -233,7 +233,7 @@
stopifnot(is.numeric(pars), length(pars) == np)
## Assumption: all parameters, including the residual SD on SD-scale
sigma <- pars[np]
- .Call(lme4Eigen:::lmer_Deviance, pp$ptr(), resp$ptr(), pars[-np]/sigma)
+ .Call(lmer_Deviance, pp$ptr(), resp$ptr(), pars[-np]/sigma)
sigsq <- sigma^2
pp$ldL2() + (resp$wrss() + pp$sqrL(1))/sigsq + n * log(2 * pi * sigsq)
}
Modified: pkg/lme4Eigen/man/mcmcsamp.Rd
===================================================================
--- pkg/lme4Eigen/man/mcmcsamp.Rd 2012-02-28 23:26:26 UTC (rev 1629)
+++ pkg/lme4Eigen/man/mcmcsamp.Rd 2012-02-29 11:57:05 UTC (rev 1630)
@@ -12,7 +12,7 @@
\item{verbose}{should verbose output be given?}
- \item{...}{}
+ \item{...}{unused currently}
}
\value{
a Markov chain Monte Carlo sample as a matrix
Modified: pkg/lme4Eigen/tests/getME.R
===================================================================
--- pkg/lme4Eigen/tests/getME.R 2012-02-28 23:26:26 UTC (rev 1629)
+++ pkg/lme4Eigen/tests/getME.R 2012-02-29 11:57:05 UTC (rev 1630)
@@ -17,7 +17,7 @@
}
fm1 <- lmer(diameter ~ (1|plate) + (1|sample), Penicillin)
-chkMod(fm1)
+chkIMod(fm1)
fm2 <- lmer(angle ~ recipe * temperature + (1|recipe:replicate), cake)
stopifnot(fixef(fm2) == getME(fm2,"beta"))
Modified: pkg/lme4Eigen/tests/nbinom.R
===================================================================
--- pkg/lme4Eigen/tests/nbinom.R 2012-02-28 23:26:26 UTC (rev 1629)
+++ pkg/lme4Eigen/tests/nbinom.R 2012-02-29 11:57:05 UTC (rev 1630)
@@ -1,5 +1,8 @@
library(lme4Eigen)
+getNBdisp <- lme4Eigen:::getNBdisp
+refitNB <- lme4Eigen:::refitNB
+
simfun <- function(sd.u=1,NBtheta=0.5,
nblock=25,
fform=~x,
@@ -15,13 +18,18 @@
}
set.seed(102)
-d1 <- simfun()
-t1 <- system.time(g1 <- glmer.nb(z~x+(1|f),data=d1,debug=TRUE))
+d.1 <- simfun()
+t1 <- system.time(g1 <- glmer.nb(z ~ x + (1|f), data=d.1, verbose=TRUE))
+g1
+## ^^ FIXME: the formula ..
-lme4Eigen:::getNBdisp(g1)
-g1B <- lme4Eigen:::refitNB(g1,theta=lme4Eigen:::getNBdisp(g1))
-deviance(g1)-deviance(g1B)
-(fixef(g1)-fixef(g1B))/fixef(g1)
+d1 <- getNBdisp(g1)
+(g1B <- refitNB(g1,theta=getNBdisp(g1)))
+(ddev <- deviance(g1)-deviance(g1B))
+(rel.d <- (fixef(g1)-fixef(g1B))/fixef(g1))
+stopifnot(all.equal(d1, 0.448, tol = 1e-3),
+ all.equal(ddev, 0.0007, tol=.0005),
+ abs(rel.d) < 0.0004)
## library(glmmADMB)
## t2 <- system.time(g2 <- glmmadmb(z~x+(1|f),
@@ -30,48 +38,64 @@
## NLL=-logLik(g2),
## theta=g2$alpha)
## 0.4487
+glmmADMB_vals <-
+ list(fixef = c("(Intercept)"=0.92871, x=2.0507),
+ NLL = structure(2944.62, class = "logLik", df= 4, nobs= 1000L),
+ theta = 0.4487)
+##' 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
+}
-glmmADMB_vals <- structure(list(fixef =
- structure(c(0.92871, 2.0507), .Names = c("(Intercept)",
-"x")), structure(2944.62, class = "logLik", df = 4, nobs = 1000L),
- 0.4487), .Names = c("fixef", "NLL", "theta"))
+stopifnot(
+ all.equal( d1, glmmADMB_vals$ theta, tol=0.0016)
+ ,
+ all.equal(fixef(g1B), glmmADMB_vals$ fixef, tol=0.01)# not so close
+ ,
+ if(FALSE) { ## df = 3 vs df = 4 --- fails! --- FIXME ??
+ all.equal(logLik.m(g1B), -glmmADMB_vals$ NLL, tol=0.001)
+ } else
+ all.equal(as.numeric(logLik.m(g1B)), as.numeric(-glmmADMB_vals$ NLL), tol= 4e-5)
+ )
-stopifnot(all.equal(glmmADMB_vals$theta,lme4Eigen:::getNBdisp(g1),
- tol=0.0016))
+if(FALSE) { ## simulation study --------------------
-## simulation study
-## 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(glmmADMB)
+ 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(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")
+ 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???)
+ 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,
@@ -85,19 +109,28 @@
## theta=g3$alpha)
glmmADMB_epil_vals <-
- structure(list(fixef = structure(c(-1.33, 0.88392, -0.92997,
- 0.47514, -0.27016, 0.33724), .Names = c("(Intercept)", "Base",
- "trtprogabide", "Age", "Visit", "Base:trtprogabide")),
- NLL = structure(624.551, class = "logLik", df = 9, nobs = 236L),
- theta = 7.4702), .Names = c("fixef", "NLL", "theta"))
+ list(fixef =
+ c("(Intercept)"= -1.33, "Base"=0.88392, "trtprogabide"=-0.92997,
+ "Age"=0.47514, "Visit"=-0.27016, "Base:trtprogabide"=0.33724),
+ NLL = structure(624.551, class = "logLik", df = 9, nobs = 236L),
+ theta = 7.4702)
-if (FALSE) {
- ## 49 seconds: too slow for regular testing!
- t4 <- system.time(g4 <- glmer.nb(y~Base*trt+Age+Visit+(Visit|subject),
- data=epil2, debug=TRUE))
- stopifnot(all.equal(glmmADMB_epil_vals$theta,lme4Eigen:::getNBdisp(g4),
- tol=0.0016))
-}
+if(!(Sys.getenv("USER") %in% c("maechler")))
+ q("no")
+## "too slow" for regular testing -- 49 (MM at lynne: 33) seconds
+(t4 <- system.time(g4 <- glmer.nb(y~ Base*trt + Age + Visit + (Visit|subject),
+ data=epil2, verbose=TRUE)))
+g4
+(Lg4 <- logLik(g4))
+attributes(Lg4) <- attributes(Lg4)[c("class","df","nobs")]
+stopifnot(
+ all.equal(getNBdisp(g4), glmmADMB_epil_vals$ theta, tol = 0.002)
+ ,
+ all.equal(fixef (g4), glmmADMB_epil_vals$ fixef, tol = 0.004)
+ ,
+ all.equal(logLik.m (g4), - glmmADMB_epil_vals$ NLL, tol = 0.0002)
+ )
+cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons''
More information about the Lme4-commits
mailing list