[Lme4-commits] r1755 - in pkg/lme4: tests testsx
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri May 25 23:43:12 CEST 2012
Author: bbolker
Date: 2012-05-25 23:43:11 +0200 (Fri, 25 May 2012)
New Revision: 1755
Added:
pkg/lme4/tests/README
Modified:
pkg/lme4/tests/nbinom.R
pkg/lme4/tests/prLogistic.R
pkg/lme4/tests/refit.R
pkg/lme4/testsx/testcolonizer.R
pkg/lme4/testsx/testcrab.R
Log:
uncomment tests where fixed PIRLS step-halving allows
more detail, analysis etc. of tests
Added: pkg/lme4/tests/README
===================================================================
--- pkg/lme4/tests/README (rev 0)
+++ pkg/lme4/tests/README 2012-05-25 21:43:11 UTC (rev 1755)
@@ -0,0 +1,13 @@
+Catalog of currently-failing examples (commented out, testsx, etc.):
+
+glmmExt.R: "fail for MM" on Gaussian/inverse examples -- seems fine for me
+lmer.R: sstudy9 example. Should *not* work; is a meaningful error message possible?
+prLogistic.R: Thailand/clustered-data example from ?prLogisticDelta example in prLogistic package
+ Presumably the problem is that 100/411 random-effect levels have only zeros -- but should this mess things up?
+ glmmML and lme4.0 give nearly identical answers
+
+profile.R: fails on CBPP profiling
+
+from testsx:
+ testcolonizer: definite case where complete separation occurs, GLM does not really give a fit
+ testcrabs: ?? not sure ??
\ No newline at end of file
Modified: pkg/lme4/tests/nbinom.R
===================================================================
--- pkg/lme4/tests/nbinom.R 2012-05-24 20:32:55 UTC (rev 1754)
+++ pkg/lme4/tests/nbinom.R 2012-05-25 21:43:11 UTC (rev 1755)
@@ -22,7 +22,7 @@
d.1 <- simfun()
t1 <- system.time(g1 <- glmer.nb(z ~ x + (1|f), data=d.1, verbose=TRUE))
g1
-## ^^ FIXME: the formula ..
+## ^^ FIXME: the formula and data results show up as ..1, ..2 ; eval.parent() etc.?
d1 <- getNBdisp(g1)
(g1B <- refitNB(g1,theta=getNBdisp(g1)))
Modified: pkg/lme4/tests/prLogistic.R
===================================================================
--- pkg/lme4/tests/prLogistic.R 2012-05-24 20:32:55 UTC (rev 1754)
+++ pkg/lme4/tests/prLogistic.R 2012-05-25 21:43:11 UTC (rev 1755)
@@ -6,3 +6,50 @@
try(glmer(rgi~ sex + pped + (1|schoolid),
data = dataset, family=binomial))
+
+if (FALSE) {
+ library(ggplot2)
+ dataset$schoolid <- factor(dataset$schoolid)
+ length(levels(dataset$schoolid))
+ ttab <- with(dataset,table(schoolid,sex,rgi,pped))
+ ttab2 <- with(dataset,table(tapply(rgi,schoolid,
+ function(x) ifelse(all(x==0),0,
+ ifelse(all(x==1),1,0.5)))))
+ ggplot(dataset,aes(x=sex,y=rgi,colour=factor(pped)))+stat_sum(alpha=0.5)
+
+ ## library(glmmML)
+ ## glmmML_fit <- glmmML(rgi~ sex + pped , cluster=schoolid,
+ ## data = dataset, family=binomial)
+ ## glmmML_est <- list(sigma=glmmML_fit$sigma,beta=coef(glmmML_fit))
+ ## dput(glmmML_est)
+ glmmML_est <- structure(list(sigma = 1.25365353546143,
+ beta = structure(c(-2.19478801858317,
+ 0.548884468743364, -0.623835613907385), .Names = c("(Intercept)",
+ "sex", "pped"))),
+ .Names = c("sigma", "beta"))
+
+ ## library(lme4.0)
+ ## lme4.0_fit <- glmer(rgi~ sex + pped + (1|schoolid),
+ ## data = dataset, family=binomial)
+ ## lme4.0_est <- list(sigma=unname(sqrt(unlist(VarCorr(lme4.0_fit)))),
+ ## beta=fixef(lme4.0_fit))
+ ## dput(lme4.0_est)
+ ## detach("package:lme4.0")
+ lme4.0_est <- structure(list(sigma = 1.25369539060849, beta = structure(c(-2.19474529099587,
+ 0.548900267825802, -0.623934772981894), .Names = c("(Intercept)",
+ "sex", "pped"))), .Names = c("sigma", "beta"))
+
+ devfun <- glmer(rgi~ sex + pped + (1|schoolid),
+ data = dataset, family=binomial,
+ devFunOnly=TRUE)
+ with(glmmML_est,devfun(c(sigma,beta)))
+ with(glmmML_est,devfun(c(sigma,beta))) ## FAILS
+ devfun <- glmer(rgi~ sex + pped + (1|schoolid),
+ data = dataset, family=binomial,
+ devFunOnly=TRUE)
+ with(lme4.0_est,devfun(c(sigma,beta))) ## 6326.456
+ g3 <- with(glmmML_est,glmer(rgi~ sex + pped + (1|schoolid),
+ data = dataset, family=binomial,
+ start=list(theta=sigma,fixef=beta)),
+ verbose=10)
+}
Modified: pkg/lme4/tests/refit.R
===================================================================
--- pkg/lme4/tests/refit.R 2012-05-24 20:32:55 UTC (rev 1754)
+++ pkg/lme4/tests/refit.R 2012-05-25 21:43:11 UTC (rev 1755)
@@ -25,11 +25,10 @@
gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
cbpp, binomial)
gm1R <- refit(gm1,with(cbpp,cbind(incidence,size-incidence)))
-## FIXME: PIRLS fail
-## gm1S <- refit(gm1,simulate(gm1)[[1]])
+gm1S <- refit(gm1,simulate(gm1)[[1]])
## FIXME: testing all-zero responses
-## this gives "step factor reduced below 0.001 without reducing pwrss"
+## this gives "pwrssUpdate did not converge in 30 iterations"
## not sure if it's pathological or not
if (FALSE) {
sim1Z <- simulate(gm1)[[1]]
@@ -40,19 +39,16 @@
stopifnot(all.equal(gm1,gm1R,tol=4e-5))
getinfo(gm1)
getinfo(gm1R)
-## FIXME:: fallout from PIRLS fail
-## getinfo(gm1S)
+getinfo(gm1S)
## binomial GLMM (prob/weights)
gm2 <- glmer(incidence/size ~ period + (1 | herd), cbpp, binomial, weights=size)
gm2R <- refit(gm2,with(cbpp,incidence/size))
-## FIXME: PIRLS failure
-## gm2S <- refit(gm2,simulate(gm2)[[1]])
+gm2S <- refit(gm2,simulate(gm2)[[1]])
stopifnot(all.equal(gm2,gm2R,tol=4e-5))
getinfo(gm2)
getinfo(gm2R)
-## FIXME:: fallout from PIRLS fail
-## getinfo(gm2S)
+getinfo(gm2S)
## Bernoulli GLMM (specified as factor)
data(Contraception,package="mlmRev")
Modified: pkg/lme4/testsx/testcolonizer.R
===================================================================
--- pkg/lme4/testsx/testcolonizer.R 2012-05-24 20:32:55 UTC (rev 1754)
+++ pkg/lme4/testsx/testcolonizer.R 2012-05-25 21:43:11 UTC (rev 1755)
@@ -1,6 +1,8 @@
library(lme4.0)
## Emacs M-<Enter> --> setwd() correctly
(load("colonizer_rand.RData"))
+m0.0 <- glm(colonizers~Treatment*homespecies*respspecies, data=randdat, family=poisson)
+with(randdat,tapply(colonizers,list(Treatment,homespecies,respspecies),sum))
summary(m1.0 <- glmer(form1, data=randdat, family=poisson))
summary(m2.0 <- glmer(form2, data=randdat, family=poisson))
@@ -12,4 +14,4 @@
try(m1 <- glmer(form1,data=randdat, family=poisson, verbose=10L))
try(m1 <- glmer(form1,data=randdat, family=poisson, verbose=10L, tolPwrss=1e-13))
try(m2 <- glmer(form2,data=randdat, family=poisson, verbose=10L))
-
+detach("package:lme4", unload=TRUE)
Modified: pkg/lme4/testsx/testcrab.R
===================================================================
--- pkg/lme4/testsx/testcrab.R 2012-05-24 20:32:55 UTC (rev 1754)
+++ pkg/lme4/testsx/testcrab.R 2012-05-25 21:43:11 UTC (rev 1755)
@@ -22,6 +22,7 @@
if (packageVersion("lme4")>"0.999375-42") {
## using development lme4 ...
try(glmer1 <- glmer(fr2,weights=initial.snail.density,family ="binomial", data=randdata))
+ ## pwrssUpdate did not converge
try(glmer1B <- glmer(fr,family ="binomial", data=randdata))
if (require("lme4.0")) {
detach("package:lme4")
@@ -29,7 +30,13 @@
glmer1 <- glmer(fr2,weights=initial.snail.density,family ="binomial", data=randdata)
## alive/dead formulation
glmer1B <- glmer(fr,family ="binomial", data=randdata)
+ coef(glmer1B)
+ detach("package:lme4.0")
}
+ if (require("glmmADMB")) {
+ ## prop/weights formulation
+ glmer1B <- glmmadmb(fr,family ="binomial", data=randdata)
+ }
} else {
## CRAN version of lme4
glmer1 <- glmer(fr2,weights=initial.snail.density,family ="binomial", data=randdata)
More information about the Lme4-commits
mailing list