Bug in afex, pbkrtest or lme4 in parametric bootstrapping?
This takes care of the proximal problem but I still have a problem
with this example with the latest versions of everything --
data_alates=read.csv("alates_raw_FINAL3.csv", header=TRUE, fill=TRUE)
## devel pbkrtest from:
## http://people.math.aau.dk/~sorenh/software/pbkrtest/devel/
## devel afex from r-forge
## devel lme4 from github
library(afex)
sapply(c("afex","lme4","pbkrtest"),function(x)
paste(packageVersion(x),collapse="."))
## afex lme4 pbkrtest
## "0.6.77" "1.1.0" "0.3.5.1"
mixed(ALATES~ANT_TENDED*MELEZITOSE+NR_START+
(1|CLONES)+(1|PLANT)+(1|REPLICATE),type=3,
family=poisson,method="PB",args.test = list(nsim = 10),
data=data_alates)
## Fitting 6 lmer() models:
## [......]
## Obtaining 5 p-values:
## [Error in pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac,
verbose) :
## Downdated VtV is not positive definite
## In addition: Warning message:
## In mixed(ALATES ~ ANT_TENDED * MELEZITOSE + NR_START + (1 | CLONES) + :
## Numerical variables NOT centered on 0 (i.e., likely bogus results
if in interactions): NR_START
It seems to fail on the first submodel, which is the model without an
intercept (a pretty weird model, if you ask me ...)
mm <- model.matrix(mm1)
form <- reformulate(c("-1",
colnames(mm)[-1],"(1|CLONES)","(1|PLANT)","(1|REPLICATE)"),response="ALATES")
newdat <- data.frame(list(model.matrix(mm1)[,-1],
getME(mm1,"flist"),
model.frame(mm1)["ALATES"]),check.names=FALSE)
mm2 <- update(mm1,formula=form,data=newdat)
confint(mm2,method="boot",nsim=20)
## 14/20 bootstrap runs failed
I don't know whether it is better to devote effort to fixing this on
the lme4 level (i.e. making the optimization more robust), or the
PBmodcomp level (i.e. making it handle fit failures better)
On 13-08-21 02:34 PM, S?ren H?jsgaard wrote:
Please use the lme4 version on github; that should fix the problem, otherwise please report back to me.
Best regards
S?ren
-----Original Message-----
From: Tom Wenseleers [mailto:Tom.Wenseleers at bio.kuleuven.be]
Sent: 21. august 2013 20:20
To: Tom Wenseleers; sorenh at mail.dk; henrik.singmann at psychologie.uni-freiburg.de; Ben Bolker (bbolker at gmail.com); r-sig-mixed-models at r-project.org
Subject: RE: [R-sig-ME] Bug in afex, pbkrtest or lme4 in parametric bootstrapping?
Oh yes, and FYI the traceback() was
10: stop(gettextf("unable to find an inherited method for function %s for signature %s",
sQuote(fdef at generic), sQuote(cnames)), domain = NA)
9: (function (classes, fdef, mtable)
{
methods <- .findInheritedMethods(classes, fdef, mtable)
if (length(methods) == 1L)
return(methods[[1L]])
else if (length(methods) == 0L) {
cnames <- paste0("\"", sapply(classes, as.character),
"\"", collapse = ", ")
stop(gettextf("unable to find an inherited method for function %s for signature %s",
sQuote(fdef at generic), sQuote(cnames)), domain = NA)
}
else stop("Internal error in finding inherited methods; didn't return a unique method",
domain = NA)
})(list("glmerMod"), function (object, newdata, ...)
standardGeneric("refit"), <environment>)
8: refit(sm, newresp = yyy)
7: .getref(largeModel, smallModel, nsim = nsim, seed = seed)
6: PBrefdist.merMod(largeModel, smallModel, nsim = nsim, cl = cl,
details = details)
5: PBrefdist(largeModel, smallModel, nsim = nsim, cl = cl, details = details)
4: PBmodcomp.merMod(largeModel = <S4 object of class "glmerMod">,
smallModel = <S4 object of class "glmerMod">, nsim = 10)
3: (function (largeModel, smallModel, nsim = 1000, ref = NULL, cl = NULL,
details = 0)
{
UseMethod("PBmodcomp")
})(largeModel = <S4 object of class "glmerMod">, smallModel = <S4 object of class "glmerMod">,
nsim = 10)
2: do.call(PBmodcomp, args = c(largeModel = full.model, smallModel = fits[[c]],
args.test))
1: mixed(use ~ age + I(age^2) + urban + livch + (1 | district),
family = binomial, data = Contraception, args.test = list(nsim = 10),
method = "PB")
So I gather the error might have to do something with the lme4 function not currently working with glmerMod classes perhaps? Is that right?
Cheers,
Tom
-----Original Message-----
From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Tom Wenseleers
Sent: 21 August 2013 19:47
To: sorenh at mail.dk; henrik.singmann at psychologie.uni-freiburg.de; Ben Bolker (bbolker at gmail.com); r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] Bug in afex, pbkrtest or lme4 in parametric bootstrapping?
Dear all,
I was just trying to use parametric bootstrapping to calculate significance levels in a Poisson mixed model using afex/pbkrtest. Using the data in attachment and the code
data_alates=read.csv("alates_raw_FINAL3.csv", header=TRUE, fill=TRUE) mixed(ALATES~ANT_TENDED*MELEZITOSE+NR_START+(1|CLONES)+(1|PLANT)+(1|REPLICATE),type=3,family=poisson,method="PB",args.test = list(nsim = 10),data=data_alates)
I get the following error message though:
Fitting 6 lmer() models:
[......]
Obtaining 5 p-values:
[
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function 'refit' for signature '"glmerMod"'
In addition: Warning message:
In mixed(ALATES ~ ANT_TENDED * MELEZITOSE + NR_START + (1 | CLONES) + :
Numerical variables NOT centered on 0 (i.e., likely bogus results if in interactions): NR_START
Any thoughts perhaps what might be wrong, and what workaround I could perhaps use to make this work?
Cheers,
Tom
PS I was using the following package versions:
sessionInfo()
R version 3.0.0 (2013-04-03)
Platform: x86_64-w64-mingw32/x64 (64-bit)
locale:
[1] LC_COLLATE=English_United Kingdom.1252 LC_CTYPE=English_United Kingdom.1252 LC_MONETARY=English_United Kingdom.1252
[4] LC_NUMERIC=C LC_TIME=English_United Kingdom.1252
attached base packages:
[1] grid stats4 splines parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] lsmeans_1.10-00 plyr_1.8 effects_2.2-4 colorspace_1.2-2 nlme_3.1-109 afex_0.6-77 stringr_0.6.2 reshape2_1.2.2
[9] car_2.0-16 nnet_7.3-6 coin_1.0-22 modeltools_0.2-19 mvtnorm_0.9-9994 survival_2.37-4 pbkrtest_0.3-5.1 lme4_1.1-0
[17] Matrix_1.0-12 lattice_0.20-15 glmmADMB_0.7.7 R2admb_0.7.5.3 MASS_7.3-26 devtools_1.2
loaded via a namespace (and not attached):
[1] digest_0.6.3 evaluate_0.4.3 httr_0.2 memoise_0.1 minqa_1.2.1 multcomp_1.2-18 RCurl_1.95-4.1 tools_3.0.0 whisker_0.3-2