Bug in afex, pbkrtest or lme4 in parametric bootstrapping?

Although I also see myself out of the discussion on how to solve the 
bigger issue, I want to make one point clear, the problem is not 
related with the somewhat controversial way of how afex construes the 
model. It rather seems to be a problem of model and/or data.

By using method = "LRT" (which is only available in the development 
version from r-forge) we can get only the models and see that 
PBmodcomp fails for all models, also the unproblematic ones (the 
p-values are then taken from anova(full.model, lower.model) which does not fail):

require(afex)
data_alates <- read.csv("alates_raw_FINAL3.csv", header=TRUE, 
fill=TRUE)

t1 <-
mixed(ALATES~ANT_TENDED*MELEZITOSE+NR_START+(1|CLONES)+(1|PLANT)+(1|RE
PLICATE), family=poisson, method="LRT", data=data_alates) ## Fitting 6 
lmer() models:
## [......]
## 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


# full model versus model without NR_START (has only main effect) 
PBmodcomp(t1$full.model, t1$restricted.models[[4]]) ## Error: 
pwrssUpdate did not converge in 30 iterations ## In addition: Warning 
messages:
## 1: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, verbose) :
##   Cholmod warning 'not positive definite' at
file:../Cholesky/t_cholmod_rowfac.c, line 431 ## 2: In pwrssUpdate(pp, 
resp, tolPwrss, GQmat, compDev, fac, verbose) :
##   Cholmod warning 'not positive definite' at
file:../Cholesky/t_cholmod_rowfac.c, line 431


# full model versus model without interaction PBmodcomp(t1$full.model, 
t1$restricted.models[[5]]) ## Error in pwrssUpdate(pp, resp, tolPwrss, 
GQmat, compDev, fac, verbose) :
##   Downdated VtV is not positive definite


Note furthermore that the exact error message varies. The "Downdated 
VtV" can be replaced with the "did not converge in 30 iterations" for 
the second comparison and the extra warnings appear from time to time.
However, even when running multiple times, I don't obtain a working 
solution (although I haven't put it in a while loop with tryCatch).

Cheers,
Henrik

Ben Bolker schrieb:

   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