[R-sig-ME] feedback: Anova (type III-tests) table based on LRT for glmmTMB models (drop1, anova, mixed)
Henrik Singmann
@ingm@nn @ending from p@ychologie@uzh@ch
Thu Aug 9 15:00:22 CEST 2018
Hi Guillaume,
Thanks for your feedback. Indeed, the initial version could not handle
fancy formulas. I have just updated the package on github to allow this
and also ensured it works with glmmTMB:
library("monet")
library("glmmTMB")
## adapted from the vignette:
set_sum_contrasts()
Owls <- transform(Owls,
Nest=reorder(Nest,NegPerChick),
NCalls=SiblingNegotiation,
FT=FoodTreatment)
zipp_test <- test_terms(formula = NCalls~(FT+scale(ArrivalTime))*SexParent +
offset(log(BroodSize)),
data = Owls, extra_formula = ~ (1|Nest),
est_fun = glmmTMB,
arg_est = list(ziformula=~1, family=poisson)
)
zipp_test
# glmmTMB Anova Table (Type III tests)
#
# Model: NCalls ~ (FT + scale(ArrivalTime)) * SexParent +
offset(log(BroodSize)) +
# Model: (1 | Nest)
# Data: Owls
# Effect Df 1 Df 0 Chisq Pr(>Chisq)
# 1 FT 8 1 20.03 *** <.0001
# 2 scale(ArrivalTime) 8 1 58.62 *** <.0001
# 3 SexParent 8 1 0.00 >.99
# 4 FT:SexParent 8 1 0.00 >.99
# 5 scale(ArrivalTime):SexParent 8 1 0.00 >.99
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1
If you have some example data, where the LRT fails, I could take a look
at it as well. Ideally, I would then add this as a test to the package.
So preferably with some data that is available in some package or so.
However, if this is not possible, I could also only run the tests only
locally and keep your data private. In this case, send me the data offlist.
Cheers,
Henrik
Am 09.08.2018 um 10:20 schrieb Guillaume Adeux:
> Hi Henrik,
>
> Thank you very much for this work. This is highly appreciated. Here is
> some feedback...
>
> I did this on my data with the following code:
>
> #set_sum_contrasts()
> #test_terms(Ratio~block+year_s+syst+syst:year_s,extra_formula=~(1|year:plot)+(1|plot)+(1|year),est_fun=glmmTMB::glmmTMB,data=bio)
> Here is the output (year_s is scale(year)):
> glmmTMB::glmmTMB Anova Table (Type III tests) Model: Ratio2 ~ block +
> year_s + syst * year_s + (1 | year:plot) + (1 | Model: plot) + (1 |
> year) Data: bio Effect Df 1 Df 0 Chisq Pr(>Chisq) 1 block 15 1 NA <NA>
> 2 year_s 15 1 0.00 .99 3 syst 15 4 10.75 * .03 4 year_s:syst 15 4 0.01
> >.99 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1
>
> -The function does not seem to recognize scale() when it is directly
> applied in the formula (i.e.
> Ratio2~block+scale(year)+syst+syst:scale(year)). The variable needs to
> be scaled outside of the formula. Otherwise, we get an error message
> "Error in model.matrix.default(formula,data=new_data): model frame and
> formula mismatch in model.matrix()". The same thing happens if you
> want to add a small constant to the response (i.e. Ratio + 0.00001).
> This needs to be done outside of the function.
>
> -When a submodel has a convergence problem, the function is not
> capable of doing the likelihood ratio test (cf. submodel in which
> "block" was removed). This is not so surprising but it seems like like
> afex::mixed was capable of doing so. I tried bumping up the number of
> iterations in arg_est with arg_est=list(control=glmmTMBControl(...)) -
> which works - but did not resolve the problem.
>
> I can send you the data if you want to have a look see :) I'm afraid I
> can't add much more.
>
> Thanks again! You have a pending citation!
>
> Guillaume ADEUX
>
>
> 2018-08-07 17:42 GMT+02:00 Henrik Singmann
> <singmann using psychologie.uzh.ch <mailto:singmann using psychologie.uzh.ch>>:
>
> Hi all,
>
> It took me some time, but I managed to come up with something that
> might be of help here. Specifically, I have worked on a new
> package, monet, that runs Type III like tests with arbitrary
> estimation function using anova() (i.e., LRT) as default. see:
> https://github.com/singmann/monet <https://github.com/singmann/monet>
> It basically generalizes what afex::mixed does for arbitrary
> estimation and test functions.
>
> I do not have glmmTMB installed, but as of now it works with lm
> and lme4::lmer. The main function is test_term(). It allows to
> pass two formulas, one for which submodels are created and
> estimates (i.e., the fixed-effects part) and one additional part
> which can for example hold the random-effects part.
>
> devtools::install_github("singmann/monet")
> library("monet")
> set_sum_contrasts() ## quite important, currently coding is not
> checked
> data("Machines", package = "MEMSS")
>
> # ignoring repeated-measures
> m1 <- test_terms(score ~ Machine, data=Machines, est_fun = lm)
> m1
> # lm Anova Table (Type III tests)
> #
> # Model: score ~ Machine
> # Data: Machines
> # Effect Df 1 Df 0 F Pr(>F)
> # 1 Machine 51 2 26.30 *** <.0001
> # ---
> # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1
>
> # simple model with random-slopes for repeated-measures factor
> m3 <- test_terms(score ~ Machine, data=Machines,
> extra_formula = ~ (Machine|Worker),
> est_fun = lme4::lmer, arg_est = list(REML = FALSE),
> arg_test = list(model.names=c("f", "r")))
> m3
> # lme4::lmer Anova Table (Type III tests)
> #
> # Model: score ~ Machine + (Machine | Worker)
> # Data: Machines
> # Effect Df 1 Df 0 Chisq Pr(>Chisq)
> # 1 Machine 10 2 17.14 *** .0002
> # ---
> # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1
>
> It returns an object of class "monet" with print(), nice(), and
> anova() methods. The package is also quite lightweight and has no
> strong dependencies besides to stats.
>
> It it does not work with glmmTMB, a reproducible example would be
> great. And I am of course happy to hear of any other comments.
>
> Cheers,
> Henrik
>
>
>
> Am 02.08.2018 um 09:22 schrieb Guillaume Adeux:
>
> Sorry, my intent was definetely not to shortcut anyone or
> anything.
>
> I thought I would give you all a little feedback on your
> propositions.
>
> The problem with drop1 (when an interaction is in the model)
> is when you
> want it to behave like a type III anova. This results in one
> of the
> variables having zero degree of freedom, hence LRT of 0 and no
> p-value.
>
> Here is the example on an output:
>
>
>
> Model: mod=glmmTMB(ratio ~ block + trt * time
> +(1|plot)+(1|year)+(1|year:plot),family=list(family="beta",link="logit"),data=biomass)
>
> drop1(mod,.~block+trt*time,test=”Chisq”)
>
> Single term deletions
>
> Df AIC
> LRT Pr(>Chi)
>
> <none> -5092.2
>
> block 1 -5093.5
> 0.7345 0.391415
>
> trt 4 -5082.7
> 17.5018 0.001544
>
> time 0 -5092.2
> 0
>
> trt:time 4 -5100.0
> 0.2169 0.994527
>
>
> If I understand correctly, drop1 is incapable of comparing a
> model A+A:B
> with A+B+A:B.
>
> Anova.III.glmmTMB indeed works but only yields Wald tests
> which I thought
> were not ideal for glmms.
>
> I contacted the developper of the {afex} package to see if the
> mixed
> function could be adapted to run on glmmTMB objects but no
> luck as of now.
> Considering the framework of glmmTMB is similar to glmer,
> maybe this can be
> easily adapted.
>
> Thanks for your interest. Don't hesitate if you have any input.
>
> Guillaume ADEUX
>
> [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-sig-mixed-models using r-project.org
> <mailto:R-sig-mixed-models using r-project.org> mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
>
>
--
Dr. Henrik Singmann
PostDoc
Universität Zürich, Schweiz
http://singmann.org
More information about the R-sig-mixed-models
mailing list