[R-sig-ME] Bug in afex, pbkrtest or lme4 in parametric bootstrapping?
Søren Højsgaard
sorenh at math.aau.dk
Thu Aug 22 00:20:38 CEST 2013
Well, the biggest room in the world is the room for improvement :) I'll look into it within the next few days and report back.
Cheers
Søren
-----Original Message-----
From: Ben Bolker [mailto:bbolker at gmail.com]
Sent: 22. august 2013 00:12
To: Henrik Singmann
Cc: Søren Højsgaard; Tom Wenseleers; r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Bug in afex, pbkrtest or lme4 in parametric bootstrapping?
Fair enough. I haven't checked PBmodcomp carefully -- but errors in fitting particular PB replicates are not surprising to me, especially with relatively small data sets. I don't know if PBmodcomp is designed to be robust to these types of failures, but it probably should be ...
That said, I'm sure there is also room for improvement in the robustness of lme4.
On 13-08-21 05:32 PM, Henrik Singmann wrote:
>
> 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
>>>
>>>
>
More information about the R-sig-mixed-models
mailing list