[Lme4-commits] r1620 - pkg/lme4Eigen/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 24 18:21:37 CET 2012
Author: bbolker
Date: 2012-02-24 18:21:36 +0100 (Fri, 24 Feb 2012)
New Revision: 1620
Modified:
pkg/lme4Eigen/tests/elston.R
pkg/lme4Eigen/tests/nbinom.R
Log:
cleaned up tests -- remove dependencies to glmmADMB, ggplot2,
get rid of longer-running stuff
Modified: pkg/lme4Eigen/tests/elston.R
===================================================================
--- pkg/lme4Eigen/tests/elston.R 2012-02-24 16:25:33 UTC (rev 1619)
+++ pkg/lme4Eigen/tests/elston.R 2012-02-24 17:21:36 UTC (rev 1620)
@@ -1,61 +1,75 @@
-tickdata <- read.table("Elston2001_tickdata.txt",header=TRUE,
- colClasses=c("factor","numeric","factor","numeric","factor","factor"))
+## original code for reading/aggregating:
-tickdata <- transform(tickdata,cHEIGHT=scale(HEIGHT,scale=FALSE))
+## tickdata <- read.table("Elston2001_tickdata.txt",header=TRUE,
+## colClasses=c("factor","numeric","factor","numeric","factor","factor"))
+## tickdata <- transform(tickdata,cHEIGHT=scale(HEIGHT,scale=FALSE))
+## for (i in names(tickdata)) {
+## if (is.factor(tickdata[[i]])) {
+## tickdata[[i]] <- factor(tickdata[[i]],levels=sort(as.numeric(levels(tickdata[[i]]))))
+## }
+## }
+## summary(tickdata)
+## grouseticks <- tickdata
+
+## library(reshape)
+## meantick <- rename(aggregate(TICKS~BROOD,data=tickdata,FUN=mean),
+## c(TICKS="meanTICKS"))
+## vartick <- rename(aggregate(TICKS~BROOD,data=tickdata,FUN=var),
+## c(TICKS="varTICKS"))
+## uniqtick <- unique(subset(tickdata,select=-c(INDEX,TICKS)))
+## grouseticks_agg <- Reduce(merge,list(meantick,vartick,uniqtick))
+
+## save("grouseticks","grouseticks_agg",file="grouseticks.rda")
+
+library(lme4Eigen)
+data(grouseticks)
+do.plots <- FALSE
form <- TICKS~YEAR+HEIGHT+(1|BROOD)+(1|INDEX)+(1|LOCATION)
+
## fit with lme4
-library(lme4)
-t1 <- system.time(full_mod1 <- glmer(form, family="poisson",data=tickdata))
-c1 <- c(fixef(full_mod1),unlist(VarCorr(full_mod1)), logLik=logLik(full_mod1),time=t1["elapsed"])
-allcoefs1 <- c(unlist(full_mod1 at ST),fixef(full_mod1))
-detach("package:lme4")
+## library(lme4)
+## t1 <- system.time(full_mod1 <- glmer(form, family="poisson",data=grouseticks))
+## c1 <- c(fixef(full_mod1),unlist(VarCorr(full_mod1)), logLik=logLik(full_mod1),time=t1["elapsed"])
+## allcoefs1 <- c(unlist(full_mod1 at ST),fixef(full_mod1))
+## detach("package:lme4")
-## fit with lme4Eigen
-library(lme4Eigen)
-t2 <- system.time(full_mod2 <- glmer(form, family="poisson",data=tickdata))
-c2 <- c(fixef(full_mod2),unlist(VarCorr(full_mod2)), logLik=logLik(full_mod2),time=t2["elapsed"])
+## lme4 summary results:
+t1 <- structure(c(1.288, 0.048, 1.36, 0, 0), class = "proc_time",
+ .Names = c("user.self",
+ "sys.self", "elapsed", "user.child", "sys.child"))
+c1 <- structure(c(11.3559080756861, 1.1804105508475, -0.978704335712111,
+ -0.0237607330254979, 0.293232458048324, 0.562551624933584,
+ 0.279548178949372,
+ -424.771990224991, 1.36),
+ .Names = c("(Intercept)", "YEAR96",
+ "YEAR97", "HEIGHT", "INDEX", "BROOD", "LOCATION",
+ "logLik", "time.elapsed"
+ ))
+allcoefs1 <- structure(c(0.541509425632023, 0.750034415832756,
+ 0.528723159081737,
+ 11.3559080756861, 1.1804105508475,
+ -0.978704335712111, -0.0237607330254979
+ ),
+ .Names = c("", "", "", "(Intercept)",
+ "YEAR96", "YEAR97", "HEIGHT"))
+
+t2 <- system.time(full_mod2 <- glmer(form, family="poisson",data=grouseticks))
+c2 <- c(fixef(full_mod2),unlist(VarCorr(full_mod2)),
+ logLik=logLik(full_mod2),time=t2["elapsed"])
+
allcoefs <- function(x) c(getME(x,"theta"),getME(x,"beta"))
## deviance function
-mm <- glmer(form, family="poisson",data=tickdata,devFunOnly=TRUE,tolPwrss=1e-12,verbose=4,
- compDev=FALSE)
-mm(allcoefs1)
-## works with compDev=FALSE, fails with compDev=TRUE
+mm <- glmer(form, family="poisson",data=grouseticks,devFunOnly=TRUE)
+## ,compDev=FALSE)
+mm2 <- glmer(form, family="poisson",data=grouseticks,
+ devFunOnly=TRUE,compDev=TRUE)
+stopifnot(all.equal(mm(allcoefs1),mm2(allcoefs1)))
## refit
-full_mod3 <- refit(full_mod2,tickdata$TICKS)
+full_mod3 <- refit(full_mod2,grouseticks$TICKS)
-fn <- "elston_fits.RData"
-## FIXME::: some possibility of differing results? 1780.
-## what changed ???
-if (!file.exists(fn)) {
- tvec <- seq(-13,-7,by=0.1)
- dmat <- matrix(nrow=length(tvec),ncol=9,
- dimnames=list(NULL,c("deviance","time_elapsed",
- paste("theta",1:3,sep=""),paste("beta",1:4,sep=""))))
- for (i in seq_along(tvec)) {
- tt <- system.time(gg <- glmer(form,family="poisson",data=tickdata,tolPwrss=10^tvec[i]))["elapsed"]
- cat(i,tvec[i],deviance(gg),"\n")
- dmat[i,] <- c(deviance(gg),tt,allcoefs(gg))
- }
- dmat <- data.frame(logtol=tvec,dmat)
- save("dmat",file=fn)
-} else load(fn)
-
-newdev <- apply(dmat[,-(1:3)],1,mm)
-library(ggplot2)
-library(reshape)
-qplot(logtol,value,data=melt(dmat,id.var="logtol"),geom=c("line"))+
- facet_wrap(~variable,scale="free")+
- geom_line(data=data.frame(logtol=dmat$logtol,value=newdev,variable="deviance"),
- colour="blue")+
- geom_hline(data=data.frame(variable="deviance",val=mm(allcoefs1)),
- aes(yintercept=val),colour="red")+theme_bw()+
- geom_vline(data=data.frame(variable="deviance",val=-10),
- aes(xintercept=val),colour="gray")
-
-detach("package:lme4Eigen")
cbind(c1,c2)
Modified: pkg/lme4Eigen/tests/nbinom.R
===================================================================
--- pkg/lme4Eigen/tests/nbinom.R 2012-02-24 16:25:33 UTC (rev 1619)
+++ pkg/lme4Eigen/tests/nbinom.R 2012-02-24 17:21:36 UTC (rev 1620)
@@ -23,16 +23,15 @@
deviance(g1)-deviance(g1B)
(fixef(g1)-fixef(g1B))/fixef(g1)
-if (FALSE) {
- library(glmmADMB)
- t2 <- system.time(g2 <- glmmadmb(z~x+(1|f),
- data=d1,family="nbinom"))
- glmmADMB_vals <- list(fixef=fixef(g2),
- NLL=-logLik(g2),
- theta=g2$alpha)
- ## 0.4487
-}
+## library(glmmADMB)
+## t2 <- system.time(g2 <- glmmadmb(z~x+(1|f),
+## data=d1,family="nbinom"))
+## glmmADMB_vals <- list(fixef=fixef(g2),
+## NLL=-logLik(g2),
+## theta=g2$alpha)
+## 0.4487
+
glmmADMB_vals <- structure(list(fixef =
structure(c(0.92871, 2.0507), .Names = c("(Intercept)",
"x")), structure(2944.62, class = "logLik", df = 4, nobs = 1000L),
@@ -41,51 +40,50 @@
stopifnot(all.equal(glmmADMB_vals$theta,lme4Eigen:::getNBdisp(g1),
tol=0.0016))
-if (FALSE) {
- ## 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)
- }
+## 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(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(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(m2,aes(x=glmer-glmmadmb))+geom_histogram()
+## ## glmer is slightly more biased (but maybe the MLE itself is biased???)
### epilepsy example:
-data(epil2,package="glmmADMB")
-epil2$subject <- factor(epil2$subject)
+data(epil,package="MASS")
+epil2 <- transform(epil,Visit=(period-2.5)/5,
+ Base=log(base/4),Age=log(age),
+ subject=factor(subject))
-if (FALSE) {
- t3 <- system.time(g3 <- glmmadmb(y~Base*trt+Age+Visit+(Visit|subject),
- data=epil2, family="nbinom"))
- glmmADMB_epil_vals <- list(fixef=fixef(g3),
- NLL=-logLik(g3),
- theta=g3$alpha)
-}
+## t3 <- system.time(g3 <- glmmadmb(y~Base*trt+Age+Visit+(Visit|subject),
+## data=epil2, family="nbinom"))
+## glmmADMB_epil_vals <- list(fixef=fixef(g3),
+## NLL=-logLik(g3),
+## 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",
More information about the Lme4-commits
mailing list