[R-sig-ME] feedback: Anova (type III-tests) table based on LRT for glmmTMB models (drop1, anova, mixed)
Guillaume Adeux
guill@ume@imon@@2 @ending from gm@il@com
Thu Aug 9 10:20:58 CEST 2018
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>:
> 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
> 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 mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list