[R-sig-ME] Chisq = 0 and p = 1 with nested lmer models
Geoff Brookshire
creepingbrain at gmail.com
Thu Sep 29 00:19:51 CEST 2011
Oops! Change
fm.F1.2a <- lmer(outcome ~ F2*F3 - F1:F2:F3 + (1 | subject) + (1 | item),
data = samp.dat, REML = FALSE)
to
fm.F1.2a <- lmer(outcome ~ F2*F3 + (1 | subject) + (1 | item),
data = samp.dat, REML = FALSE)
Sorry for the confusion!
Geoff
On Wed, Sep 28, 2011 at 6:13 PM, Geoff Brookshire
<broog731 at newschool.edu> wrote:
> Hi again,
>
> After an education-filled delay, this makes perfect sense. I. White
> pointed out to me that:
> "it will often make sense to compare a quadratic
>
> b0 + b1*x + b2*x^2
>
> with b0 + b1*x, but when would it make sense to compare it with
>
> b0 + b2*x^2 ?
>
> It all has to do with marginality: 'interactions should not be included
> without main effects, nor higher-degree polynomial terms without their
> lower-degree relatives' (McCullagh and Nelder, Generalized Linear Models
> 2nd ed p89). See also Nelder (1977) J Roy Stat Soc A 140 48-77,
> especially p49-50."
>
>
> This leads to a further question -- what's the recommended way to get
> at those lower-order effects? One easy way is to just use pvals.fnc
> from languageR on the full model.
>
> Alternatively, we could compare a fitted model containing on the main
> effect of interest (and random effects) against a model containing
> only an intercept, like so:
> fm.F1.1 <- lmer(outcome ~ F1 + (1 | subject) + (1 | item),
> data = samp.dat, REML = FALSE)
> fm.F1.2 <- lmer(outcome ~ 1 + (1 | subject) + (1 | item),
> data = samp.dat, REML = FALSE)
> anova(fm.F1.1, fm.F1.2)
>
> ... or compare two models additionally containing everything of
> interest that doesn't involve that main effect, like so:
> fm.F1.1a <- lmer(outcome ~ F1 + F2*F3 + (1 | subject) + (1 | item),
> data = samp.dat, REML = FALSE)
> fm.F1.2a <- lmer(outcome ~ F2*F3 - F1:F2:F3 + (1 | subject) + (1 | item),
> data = samp.dat, REML = FALSE)
> anova(fm.F1.1a, fm.F1.2a)
>
> The second of these model comparisons gives a p-value much closer to
> the pMCMC output by pvals.fnc, but both are a bit more significant.
> Does anyone have a recommendation about which of these is most
> trustworthy / statistically valid?
>
> Thanks,
> Geoff
>
> On Fri, Sep 16, 2011 at 6:18 PM, Ben Bolker <bbolker at gmail.com> wrote:
>> Geoff Brookshire <creepingbrain at ...> writes:
>>
>>>
>>> Dear all,
>>>
>>> I'm using lme4 to analyze some reaction time data I've collected, and
>>> am running into an odd problem. I'm not certain if it's a bug or a
>>> problem with my thinking, but hopefully (and probably) it's the latter.
>>>
>>> When trying to compare two models, one of which is the same as the
>>> other but with one term removed, the result is that there is no
>>> difference between the models: Chisq = 0, Chisq Df = 0, and p = 1.
>>>
>>> A little about the design:
>>> This is a reaction-time experiment with a 2x2x2 design over factors
>>> F1, F2, and F3. F1 is varied within-subjects and within-items, and F2
>>> & F3 are varied within-subjects and between-items.
>>>
>>> So far, I've been testing the models like this:
>>>
>>> fm.1 <- lmer(outcome ~ F1*F2*F3 + (1 | subject) + (1 | item),
>>> data = samp.dat, REML = FALSE)
>>> fm.2 <- lmer(outcome ~ F1*F2*F3 - F1:F2:F3 + (1 | subject) + (1 | item),
>>> data = samp.dat, REML = FALSE)
>>> anova(fm.1, fm.2)
>>>
>>> This works, and gives what we expect. When I try to test the 2-term
>>> interactions, though, the error comes up:
>>>
>>> fm.3 <- lmer(outcome ~ F1*F2*F3 - F1:F2 + (1 | subject) + (1 | item),
>>> data = samp.dat, REML = FALSE)
>>> anova(fm.1, fm.3)
>>>
>>> This gives the following output:
>>>
>>> fm.1: outcome ~ F1 * F2 * F3 + (1 | subject) + (1 | item)
>>> fm.3: outcome ~ F1 * F2 * F3 - F1:F2 + (1 | subject) + (1 | item)
>>> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
>>> fm.1 11 -677.34 -617.72 349.67
>>> fm.3 11 -677.34 -617.72 349.67 0 0 1
>>>
>>> The same thing happens when looking at F2:F3 and F1:F3. Any ideas?
>>> Just entering anova(fm.1) gives sensible results and F-values for the
>>> 2-way interactions.
>>
>> When you subtract F1:F2 from F1*F2*F3 you are still
>> going to get a fully specified model, it will just be reparameterized
>> (i.e. you won't reduce the number of parameters this way). This is
>> because of the way that R interprets model formulae (i.e. it's more
>> basic than lme4).
>>
>> Here's some an experiment confirming this:
>>
>>> d <- expand.grid(F1=LETTERS[1:2],F2=letters[1:2],F3=factor(1:2))
>>> d
>> F1 F2 F3
>> 1 A a 1
>> 2 B a 1
>> 3 A b 1
>> 4 B b 1
>> 5 A a 2
>> 6 B a 2
>> 7 A b 2
>> 8 B b 2
>>> model.matrix(~F1*F2*F3,d)
>> (Intercept) F1B F2b F32 F1B:F2b F1B:F32 F2b:F32 F1B:F2b:F32
>> 1 1 0 0 0 0 0 0 0
>> 2 1 1 0 0 0 0 0 0
>> 3 1 0 1 0 0 0 0 0
>> 4 1 1 1 0 1 0 0 0
>> 5 1 0 0 1 0 0 0 0
>> 6 1 1 0 1 0 1 0 0
>> 7 1 0 1 1 0 0 1 0
>> 8 1 1 1 1 1 1 1 1
>>
>>
>>> model.matrix(~F1*F2*F3-F1:F2,d)
>> (Intercept) F1B F2b F32 F1B:F32 F2b:F32 F1B:F2b:F31 F1B:F2b:F32
>> 1 1 0 0 0 0 0 0 0
>> 2 1 1 0 0 0 0 0 0
>> 3 1 0 1 0 0 0 0 0
>> 4 1 1 1 0 0 0 1 0
>> 5 1 0 0 1 0 0 0 0
>> 6 1 1 0 1 1 0 0 0
>> 7 1 0 1 1 0 1 0 0
>> 8 1 1 1 1 1 1 0 1
>>
>> If you have really thought this through and know that you
>> have a sensible hypothesis corresponding to setting parameters
>> 5 through 7 of the original parameterization to zero, you can do
>> it ...
>>
>> d <- transform(d,F3way=as.numeric((F1:F2:F3=="B:b:2")))
>> model.matrix(~F1+F2+F3+F3way,data=d)
>>
>> A more tedious but general way to do this is to take the
>> original model matrix, pull out the bits you want, rename
>> them if necessary (colons in names will not play nicely
>> with formulae) and use the numeric dummy variables as
>> predictors in a new model.
>>
>> Hope that helps.
>>
>> Ben Bolker
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
More information about the R-sig-mixed-models
mailing list