[R-sig-ME] degrees of freedom in Gamma GLMMs

Henrik Singmann singmann at psychologie.uzh.ch
Mon Feb 19 13:53:52 CET 2018


Dear Julia, (please keep the list in CC)

 From your mail it was not clear you reported t-tests. And in fact, the 
tests reported by glmer() per default are not really t-test but rather 
z-tests. This can be seen by the header of the p-value column. Have a 
look at an example:

summary(m1s)
# [...]
# Fixed effects:
#                 Estimate Std. Error t value Pr(>|z|)
# (Intercept)      1.05470    0.03709  28.438  < 2e-16 ***
# task1           -0.15399    0.03741  -4.116 3.86e-05 ***
# stimulus1       -0.08703    0.01325  -6.567 5.14e-11 ***
# task1:stimulus1 -0.05180    0.01385  -3.741 0.000184 ***

The p-values can be replicated via the normal distribution (pnorm):

2*pnorm(abs(fixef(m1s$full_model)) /
           sqrt(diag(vcov(m1s$full_model))), lower.tail = FALSE)
#   (Intercept)           task1       stimulus1 task1:stimulus1
# 6.837380e-178    3.858592e-05    5.141127e-11    1.835465e-04

So one option would be to tell the editors that you can use either 
z-tests or likelihood-ratio tests. In case they do not agree with this 
your best bet might be to follow Ben Bolker's advise and ask them how 
they suggest you calculate them (because it is not really clear how you 
could).

Hope that helps,
Henrik

Am 19.02.2018 um 12:20 schrieb Klaus, J. (Jana):
> Dear Henrik,
>
> Thanks a lot for your response. However, I'm not sure whether this is what the editor wants, as he requested " When you report t-tests, please add the df, using the following format: t(df) = x.xx, p = 0.xx." We responded that no consensus exists as of yet (as opposed to LMMs, even though I understand it is questionable there as well, but can at least be computed with lmerTest) and referenced the Lo & Andrews (2015) paper, to which they said " We found several citations, more recent than 2015,  to counter the Lo and Andrews (2015) position about degrees of freedom in generalized mixed models.  We think that Readers will appreciate knowing fuller details about your statistical approach.  We will also consult with statisticians about this matter, and we may correspond with you." Needless to say they haven't corresponded with me yet...
>
> Also, I'm not sure whether it is feasible to use afex for my dataset, as I only have 32 subjects and 6224 observations - from what I understand from the manual, LRT in afex are only useful with more Ns/observations, no?
>
> Thanks so much,
> Jana
>
> -----Original Message-----
> From: Henrik Singmann [mailto:singmann at psychologie.uzh.ch]
> Sent: maandag 19 februari 2018 11:27
> To: Klaus, J. (Jana) <J.Klaus at psych.ru.nl>
> Subject: !!!!!SPAM!!!!! Re: degrees of freedom in Gamma GLMMs
>
> Hi Jana,
>
> If you want degrees of freedom, you could use package afex which allows
> you to calculate likelihood-ratio tests for the fixed effects. This way
> you can simply report chi-square tests with df equal to the number of
> omitted parameters for each test. For example:
>
> library("afex")
> data("fhch2010") # load Freeman, Heathcote, Chalmers, and Hockley (2010)
> data
> fhch <- droplevels(fhch2010[ fhch2010$correct,]) # remove errors
>
> m1s <- mixed(rt ~ task*stimulus + (stimulus|id) + (task|item),
>                fhch, method = "LRT", family = Gamma(link = "identity"))
> m1s
> # Mixed Model Anova Table (Type 3 tests, LRT-method)
> #
> # Model: rt ~ task * stimulus + (stimulus | id) + (task | item)
> # Data: fhch
> # Df full model: 11
> #          Effect df     Chisq p.value
> # 1          task  1 14.43 ***   .0001
> # 2      stimulus  1 28.37 ***  <.0001
> # 3 task:stimulus  1 12.11 ***   .0005
> # ---
> # Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1
> (For simplicity I have omitted convergence warnings from the output.)
>
> Note that the values of the test statistics reported can be different
> from the ones reported by glmer directly, because they are based on
> comparing a full model against a reduced model in which each effect is
> omitted. However, models are also estimated with lme4::glmer.
>
> I would report the main effect of task as: chi^2(1) = 14.43, p = .0001
> (replace chi with the greek symbol and ^2 with power of two).
>
> Hope that helps,
> Henrik
>
>
> Am 19.02.2018 um 10:25 schrieb Klaus, J. (Jana):
>> Hi all,
>>
>>    
>>
>> I have fitted a GLMM with a Gamma distribution with crossed random effects in lme4 for reaction times, like this:
>>
>>    
>>
>> model = glmer(RT~Factor1*Factor2+ (1+Factor1+Factor2|subj) + (1+Factor1|item), data=xdat, family = Gamma(link = "identity"), control=glmerControl(optimizer = "bobyqa"))
>>
>>    
>>
>> The editor now insists on reporting degrees of freedom. From what I have found online, this is not trivial, and potentially not at all possible with the design. I fit the same model with glmmPQL but I cannot reproduce the results, presumably because the random effects are treated as nested rather than crossed there. Also, as far as I can tell, all other workarounds to compute dfs do not apply to this specific design. Nevertheless, I’m afraid I might be overlooking something essential already implemented in lme4 (or any other package for that matter). If this is not the case, could you point me towards any work that explains *why* it can’t be done for GLMMs? Any help is greatly appreciated!
>>
>>    
>>
>> Cheers,
>>
>> Jana
>>
>>
>> 	[[alternative HTML version deleted]]
>>
>
>

-- 
Dr. Henrik Singmann
PostDoc
Universität Zürich, Schweiz
http://singmann.org



More information about the R-sig-mixed-models mailing list