

   Dear Geoff, dear list:

   Sorry to recover this thread so long after, but I was writing my chapter before and now I'm going to calculate the significance of these main effects that I'm not including in my discussion anyway.

   So, I just wanted to make sure I have understood the procedure you kindly told me, but with my own data: 


   I have my minimum adequate model, which is:

   >s.hatchmodel1.103<-lmer(prophatch1~sex+inslarge+tlength+briventral+brithr+sex:brithr+tlength:brithr+briventral:brithr+(1|site/pair)+(1|year),REML=FALSE)

   and I want to test the effect of sex, tlength, briventral and brithr when they're not taking part in any interactions. The first possible method you suggest is using the pvals.fnc from languageR on the full model. I don't know wether you mean the results of the Wald's test of the summary table, but if so, that's not a possibility since lmer does not give p-values for models with a Gaussian distribution, but only for GLMMs.

   According to the second possibility, it would be simply like this (for example, to test brithr):

   >brithrmodel1<-lmer(prophatch1~brithr+(1|site/pair)+(1|year),REML=FALSE)

   >brithrmodel2<-lmer(prophatch1~1+(1|site/pair)+(1|year),REML=FALSE)

   >anova(brithrmodel1,brithrmodel2)


   And according to the third suggestion, the one I'm really interested in understanding (as you say it gives the most accurate pvals): basically, for calculating the significance of brithr, I have to remove all the interactions in which brithr is taking part (in this case, all of them), and test the model with the ramaining terms PLUS brithr against the same model without brithr, like this:

>brithrmodel1<-lmer(prophatch1~sex+inslarge+tlength+briventral+brithr+(1|site/pair)+(1|year),REML=FALSE)

>brithrmodel2<-lmer(prophatch1~sex+inslarge+tlength+briventral+(1|site/pair)+(1|year),REML=FALSE)

>anova(brithrmodel1,brithrmodel2)   


For briventral, I only have to remove one of the 2-way interactions:


>briventralmodel1<-lmer(prophatch1~sex+inslarge+tlength+briventral+brithr+sex:brithr+tlength:brithr+(1|site/pair)+(1|year),REML=FALSE)


>briventralmodel2<-lmer(prophatch1~sex+inslarge+tlength+brithr+sex:brithr+tlength:brithr+(1|site/pair)+(1|year),REML=FALSE)

>anova(briventralmodel1,briventralmodel2)

Something similar for tlength:


>tlengthmodel1<-lmer(prophatch1~sex+inslarge+tlength+briventral+brithr+sex:brithr+briventral:brithr+(1|site/pair)+(1|year),REML=FALSE)

>tlengthmodel2<-lmer(prophatch1~sex+inslarge+briventral+brithr+sex:brithr+briventral:brithr+(1|site/pair)+(1|year),REML=FALSE)


>anova(tlengthmodel1,tlengthmodel2)

And the same for sex:


>sexmodel1<-lmer(prophatch1~sex+inslarge+tlength+briventral+brithr+tlength:brithr+briventral:brithr+(1|site/pair)+(1|year),REML=FALSE)

>sexmodel2<-lmer(prophatch1~inslarge+tlength+briventral+brithr+tlength:brithr+briventral:brithr+(1|site/pair)+(1|year),REML=FALSE)


>anova(briventralmodel1,briventralmodel2)


Is that right? I hope I have understood it properly.

Thank you so much!!
 



________________________________
De: Geoff Brookshire <creepingbrain@gmail.com>
Para: r-sig-mixed-models@r-project.org
Enviado: jueves 29 de septiembre de 2011 0:19
Asunto: Re: [R-sig-ME] Chisq = 0 and p = 1 with nested lmer models

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)

[[elided Yahoo spam]]

Geoff

On Wed, Sep 28, 2011 at 6:13 PM, Geoff Brookshire
<broog731@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@gmail.com> wrote:
>> Geoff Brookshire <creepingbrain@...> 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@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>

_______________________________________________
R-sig-mixed-models@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
	[[alternative HTML version deleted]]


