Thanks for your feedback. Indeed, the initial version could not handle
fancy formulas. I have just updated the package on github to allow this
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1
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.
locally and keep your data private. In this case, send me the data offlist.
Post by Guillaume AdeuxHi Henrik,
Thank you very much for this work. This is highly appreciated. Here is
some feedback...
#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)
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
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
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
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.
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]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
<https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>