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

Mon Feb 19 17:06:05 CET 2018

If/when you find out the contradictory citations from the editor,
please forward them to the mailing list -- many of us would be
interested!  From your design, I'm guessing that the "minimal" df
would be (n-1)=31 for Factor2 and min(# items - 1, 31) for Factor1
(that assumes Factor2, Factor1 are numeric or 2-level categorical
variables).

On Mon, Feb 19, 2018 at 7:53 AM, Henrik Singmann
<singmann at psychologie.uzh.ch> wrote:
> 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
>
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models