[R-lang] Re: Main effects of categorical predictors in lmer
T. Florian Jaeger
tiflo@csli.stanford.edu
Mon Oct 10 13:07:00 PDT 2011
Hi Hossein,
so, let's take this one step at a time. Here's what you can do to answer
whether A or B or both have an effect (I am going to assume that the random
effect structure you used, (1 + A | sub) + (1 | item), has been shown to be
the appropriate one or at least sufficient; if not adapt accordingly):
l.full = lmer(response ~ A + B + A:B + (1 + A | sub) + (1 | item), data,
family="binomial")
l.AB = lmer(response ~ A + B + (1 + A | sub) + (1 | item), data,
family="binomial")
l.A = lmer(response ~ A + (1 + A | sub) + (1 | item), data,
family="binomial")
l.B = lmer(response ~ B + (1 | sub) + (1 | item), data, family="binomial")
if anova(l.full, l.A) is significant, B has an effect (main effect or
interaction).
if anova(l.full, l.B) is significant, A has an effect (main effect or
interaction).
(If there are not significant, non-significance of A or B, respectively,
does *not* follow). Given what you've shown us already, both will be
significant. Then
if anova(l.full, l.AB) is significant, you have a significant interaction of
A and B.
This will also be the case given what you've shown us so far. You don't
technically need to do this. But I wanted to let you know what you would
need to do to get the closest equivalent to the so-called omnibus test that
the ANOVA would provide.
Now, the regression model (if appropriately coded) also gives
you interpretable information about the contrasts between the three levels
of B. Which coding I would a priori use is hard to tell since you chose to
abstract away from your real data (A and B, rather than NP and .... ;)).
Btw, it usually helps if you describe your actual data and provide e;g.
head(data[,relevantcolumns]).
With regard to your email:
> contrasts(B)
> [,1] [,2]
> B2 0 1
> B1 -1 -1
> B3 1 0
>
> The three levels of B are coded with names rather than numbers and the
> baseline condition of B alphabetically comes second; that's why I manually
> recoded the matrix so that the baseline receives -1. However, as you can
> see, this matrix is different from the one Florian posted. I used the
> following code to create the contrasts for B (but I just read Ian's message
> which says this is wrong):
>
> > contrasts (B)<- contr.sum(3)
> > contrasts(B) <- cbind(c(0,-1,1),c(1,-1,0))
>
so here the second row overrides the first! Here's what I recommend:
# order of levels argument matters!
B = factor(B, levels = ("B1", "B2", "B3"))
contrasts(B) = contr.helmert(3) / 2
# dividing by 2 to give you distances of 1 between group means
which will give you:
> contrasts(B)
B1 -.5 -.5
B2 .5 -.5
B3 0 1
which is what you seem to want. Contrast 'B1' will contrast B==2 against the
baseline B==1. Contrast 'B2' will contrast B==3 against the mean of B==1 and
B==2.
I recommend you contrast code A:
A = factor(A, levels = ("baseline", "notbaseline"))
contrasts(B) = contr.sum(2) / 2
> contrasts(A)
A1 -.5
A2 .5
Then l.full (above) should be rather interpretable.
Florian
>
> So, if I want to compare the baseline condition of B with the mean of other
> two conditions of B, I need to create contrasts like Florian's? When I
> helmet code the levels of B with the baseline receiving 2 and the other two
> levels of B receiving -1, I get the following confusing contrasts. I have no
> idea what this means. I’m just hoping that I’m comparing the baseline with
> the mean of the other B levels:
> > contrasts (helmertB)
> 2
> -1 0
> 2 1
>
>
> And here is the complete output of the first model (where A is numerical
> and centered):
>
>
> Generalized linear mixed model fit by the Laplace approximation
> Formula: response ~ A * B + (A + 1 | sub) + (1 | item)
> Data: mydata
> AIC BIC logLik deviance
> 745.1 794.6 -362.6 725.1
> Random effects:
>
> Groups Name Variance Std.Dev. Corr
> item (Intercept) 0.93903 0.96903
> sub (Intercept) 2.39079 1.54622
> A 0.34163 0.58449 0.551
> Number of obs: 1038, groups: item, 42; sub, 36
>
> Fixed effects:
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) 1.1661 0.3285 3.550 0.000385 ***
> A -2.0741 0.1690 -12.270 < 2e-16 ***
> B1 -0.4166 0.1580 -2.636 0.008389 **
> B2 -0.2867 0.1783 -1.608 0.107842
> A:B1 0.9881 0.1594 6.199 5.68e-10 ***
> A:B2 -0.7101 0.1750 -4.059 4.94e-05 ***
> ---
>
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Correlation of Fixed Effects:
> (Intr) NP lngt.1 lngt.2 NP:l.1
> A 0.121
> B1 -0.173 0.251
> B2 0.003 0.060 -0.281
> A:B1 0.141 -0.304 -0.259 -0.095
> A:B2 0.021 0.056 -0.081 -0.300 -0.334
>
> Florian, you already interpreted this output but it was based on the your
> assumption of the contrasts matrix for B…
>
>
> By the way, I didn't quite get Ian's suggestion for A. Should I keep it
> numerical and centered?
>
> Many thanks,
>
> hossein
>
>
> On Mon, Oct 10, 2011 at 6:11 PM, Finlayson, Ian <IFinlayson@qmu.ac.uk>wrote:
>
>> **
>>
>> You probably are thinking in the ANOVA way (will leave that for others to
>> decide if that's good/bad, appropriate/innappropriate) but anova(model) will
>> give you an ANOVA-style table for your model (minus the p values).
>>
>> Not entirely sure what you mean about having added NP, I assume you added
>> it as another fixed effect, but it's not listed in your results, so you may
>> have meant something else. If B2 becomes significant after including NP then
>> that just suggests to me that NP cleans up noise in your data, allowing B2
>> to become significant (can't help but notice that the p value in the first
>> model is only just over .1, so it's not massively far from significance).
>>
>> I should have made the point earlier, but I will now: If you're actually
>> interested in whether there's a difference between B1 and B2, and between B1
>> and B3, then sum coding isn't appropriate. Sum coding will tell you if
>> there's a difference between the levels of the factor and the grand mean. If
>> you actually want to compare B1 with B2, and B3, then you'll need to use
>> treatment coding. Your model suggests this is experimental data (Subject,
>> Item etc.) so I'll assume it's a roughly balanced design and collinearity
>> might not be a huge concern.
>>
>> You're correct that you could use subsetting (dividing the data according
>> to A makes sense to me, but that is obviously atheoretical as I have no idea
>> what A represents) to get a clear idea of your pattern of results. Although,
>> as I've just noticed Florian showing, you can get a good picture if you
>> understand the coding and work through the coefficients.
>>
>> HTH,
>>
>> Ian
>>
>>
>>
>> -----Original Message-----
>> From: hossein karimi [mailto:karimihussein@gmail.com<karimihussein@gmail.com>
>> ]
>> Sent: Mon 10/10/2011 16:15
>> To: Finlayson, Ian
>> Cc: ling-r-lang-l@mailman.ucsd.edu
>> Subject: Re: [R-lang] Main effects of categorical predictors in lmer
>>
>> Hi Ian,
>>
>> Thank you for your help.I did what you suggested and got the following
>> output:
>> Fixed effects:
>> Estimate Std. Error z value Pr(>|z|)
>> (Intercept) 1.1661 0.3285 3.550 0.000385 ***
>> A -2.0741 0.1690 -12.270 < 2e-16 ***
>> B1 -0.4166 0.1580 -2.636 0.008389 **
>> B2 -0.2867 0.1783 -1.608 0.107842
>> A:B1 0.9881 0.1594 6.199 5.68e-10 ***
>> A:B2 -0.7101 0.1750 -4.059 4.94e-05 ***
>>
>> I think this is telling me there is a main effect of A, but my problem is
>> that I also need to report a main effect of B, and that's something I
>> don't
>> find in this output. I might be thinking in the old ANOVA way, but I think
>> I
>> should report main effects as well as a main interaction (if there is any)
>> and then if there's a significant interaction, I should proceed with
>> breaking down the main interaction to see where exactly it lies, am I
>> right?
>> or, are lmer results reported differently?
>>
>> By the way, I also ran a model with NP as a factor too and got the
>> following
>> result:
>> Fixed effects:
>> Estimate Std. Error z value Pr(>|z|)
>> (Intercept) 0.9371 0.3313 2.829 0.00468 **
>> A1 2.0857 0.1700 12.269 < 2e-16 ***
>> B1 -0.3075 0.1544 -1.991 0.04646 *
>> B2 -0.3651 0.1735 -2.105 0.03533 *
>> A1:B1 -0.9936 0.1603 -6.199 5.68e-10 ***
>> A1:B2 0.7139 0.1759 4.058 4.96e-05 ***
>>
>> The effect of B2 is significant in the first model above but not in
>> second!!! why is that?
>>
>> Many thanks,
>>
>> Hossein
>>
>>
>>
>>
>> On Mon, Oct 10, 2011 at 3:27 PM, Finlayson, Ian <IFinlayson@qmu.ac.uk
>> >wrote:
>>
>> > Hi Hossein,****
>> >
>> > ** **
>>
>> >
>> > As far as I can see, neither model is correct. If B is a predictor with
>> > three levels then the model should contain 6 fixed effects: An
>> intercept, A,
>> > B1 compared to B2, B1 compared to B3, B1 compared to B2 interacting with
>> A,
>> > and B1 compared to B3 interacting with A.****
>> >
>> > ** **
>>
>> >
>> > What I suspect you've done is actually enter B as a continuous predictor
>> > with a mean of zero (since you used scale()). You'll need to keep B as a
>> > factor, but give it appropriate contrasts:****
>> >
>> > ** **
>> >
>> > B <- factor(B) # Make sure it's a factor****
>>
>> >
>> > contrasts(B) <- contr.sum(3) # Generate a 3 level contrast matrix using
>> sum
>> > coding; alternative contr.helmert(3) for Helmhert coding****
>> >
>> > ** **
>> >
>> > Then build your models.****
>> >
>> > ** **
>>
>> >
>> > It's possible that it won't set the right level as a baseline (the
>> baseline
>> > will be whichever is a row of -1). If so then you'll need to manually
>> set up
>> > the matrix. For sum coding, something like:****
>> >
>> > ** **
>>
>> >
>> > contrasts(B) <- cbind(c(-1,0,1),c(-1,1,0)) # Assuming that B1 is the
>> first
>> > row in contrasts(B) - double check!****
>> >
>> > ** **
>>
>> >
>> > I'd also do something similar with A so that your estimates are
>> comparable
>> > (if you've centered 0 and 1, then it'll come out -.5 and .5, rather than
>> -1
>> > and 1)****
>> >
>> > ** **
>> >
>> > HTH,****
>> >
>> > ** **
>> >
>> > Ian****
>> >
>> > ** **
>> >
>> > ** **
>> >
>> > ** **
>> >
>> > ** **
>> >
>> > *From:* ling-r-lang-l-bounces@mailman.ucsd.edu [mailto:
>> > ling-r-lang-l-bounces@mailman.ucsd.edu] *On Behalf Of *hossein karimi
>> > *Sent:* 10 October 2011 15:05
>> > *To:* ling-r-lang-l@mailman.ucsd.edu
>> > *Subject:* [R-lang] Main effects of categorical predictors in lmer****
>> >
>> > ** **
>> >
>> > Dear R users,****
>>
>> >
>> > I'm using mixed effects models (lmer) to predict a binary dependent
>> > variable as a function of 1.a categorical predictor (A)with 2 levels (A1
>> and
>> > A2) , 2. another categorical predictor (B) with three levels (B1, B2 and
>> B3)
>> > and 3. The interaction between these two predictors. I have tried two
>> models
>> > but they return different results and I'm not sure which one is correct.
>> I'm
>> > interested in the main effect of B and the interaction between A and B
>> > (because A alone has a significant effect in both models). My problem is
>> > that there seem to be two sensible ways of examining the main effect of
>> B:
>> > 1. to helmert code and 2. to center. But these two methods produce
>> opposite
>> > results! I don't know which one I should use. Here are the two models
>> with
>> > some details and their outputs:****
>> >
>> > ** **
>>
>> >
>> > Model 1: 'A' is centered. 'B' is helmert coded ('B1'(baseline)=2,
>> 'B2'=-1,
>> > 'B3'=-1) so that I can get a main effect of B by checking to see whether
>> > baseline condition in B differs from the mean of B1 and B2 . The lmer
>> output
>> > returns a significant effect of B and no significant AxB interaction.
>> > However, as is highlighted below (in pink), the correlation between B
>> and
>> > the 'AxB' interaction is high (-54%). ****
>> >
>> > ****
>>
>> >
>> > > model.1<-lmer(response~A*B+(A+1|sub)+(1|item), mydata,
>> family="binomial")
>> > ****
>> >
>> > > print(model.1)****
>> >
>> > Generalized linear mixed model fit by the Laplace approximation ****
>> >
>> > Formula: response ~ A * B+ (A + 1 | sub) + (1 | item) ****
>> >
>> > Data: mydata****
>> >
>> > AIC BIC logLik deviance****
>> >
>> > 783 822.6 -383.5 767****
>> >
>> > Random effects:****
>> >
>> > Groups Name Variance Std.Dev. Corr ****
>> >
>> > item (Intercept) 0.7293 0.85399 ****
>> >
>> > sub (Intercept) 2.0871 1.44468 ****
>> >
>> > A 1.3812 1.17524 0.562 ****
>> >
>> > Number of obs: 1038, groups: item, 42; sub, 36****
>> >
>> > ****
>> >
>> > Fixed effects:****
>> >
>> > Estimate Std. Error z value Pr(>|z|) ****
>> >
>> > (Intercept) 1.05261 0.30283 3.476 0.000509 *******
>> >
>> > A -3.91080 0.32239 -12.131 < 2e-16 *******
>> >
>> > B 0.36128 0.09751 3.705 0.000211 *******
>> >
>> > A:B -0.29638 0.18681 -1.586 0.112626 ****
>> >
>> > ---****
>> >
>> > Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ****
>> >
>> > ****
>> >
>> > Correlation of Fixed Effects:****
>> >
>> > (Intr) A B****
>> >
>> > A 0.155 ****
>> >
>> > B 0.160 -0.278 ****
>> >
>> > A:B -0.156 0.238 -0.540****
>> >
>> > ****
>>
>> >
>> > Model 2: 'A' and 'B' are both centered. The lmer output returns no
>> > significant effect of B but the A:B interaction is significant. The
>> > correlations between predictors are generally lower and the correlation
>> > between B and A:B is reduced to -26%. ****
>> >
>> > ** **
>> >
>> > Generalized linear mixed model fit by the Laplace approximation ****
>> >
>> > Formula: resonse ~ A * B + (A + 1 | sub) + (1 | item) ****
>> >
>> > Data: mydata ****
>> >
>> > AIC BIC logLik deviance****
>> >
>> > 756.1 795.7 -370.1 740.1****
>> >
>> > Random effects:****
>> >
>> > Groups Name Variance Std.Dev. Corr ****
>> >
>> > item (Intercept) 0.87028 0.93289 ****
>> >
>> > sub (Intercept) 2.41707 1.55469 ****
>> >
>> > A 1.23669 1.11206 0.533 ****
>> >
>> > Number of obs: 1038, groups: item, 42; sub, 36****
>> >
>> > ****
>> >
>> > Fixed effects:****
>> >
>> > Estimate Std. Error z value Pr(>|z|) ****
>> >
>> > (Intercept) 1.1004 0.3239 3.398 0.000679 *******
>> >
>> > A -4.0941 0.3248 -12.605 < 2e-16 *******
>> >
>> > B -0.1461 0.1400 -1.043 0.296851 ****
>> >
>> > A:B 1.7923 0.2818 6.360 2.01e-10 *******
>> >
>> > ---****
>> >
>> > Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ****
>> >
>> > ****
>> >
>> > Correlation of Fixed Effects:****
>> >
>> > (Intr) A B****
>> >
>> > A 0.138 ****
>> >
>> > B -0.148 0.185 ****
>> >
>> > A:B 0.106 -0.292 -0.265****
>> >
>> > ****
>>
>> >
>> > I personally think Model 2 is better but the thing is that I have
>> centered
>> > a categorical predictor with *three* levels. In my searches in the web,
>> I
>>
>> > have never seen a three-level predictor to be centered; they were all
>> > two-level categorical predictors. ****
>>
>> >
>> > I have used the scale() function to center the predictors (I first
>> > converted them to numeric variables and then used the scale () function
>> to
>> > center them). As I mentioned, my problem is that I don't know how to get
>> a
>> > main effect of B as well as a *main* A:B interaction. On the one hand,
>> it
>>
>> > seems logical to compare 'B1' (baseline) with the mean of the other two
>> B
>> > conditions to see if the B manipulation has a general effect. On the
>> other
>> > hand, I hear that one needs to center variables to get a main effect.
>> ****
>> >
>> > ** **
>> >
>> > I would be grateful of you could please help me****
>> >
>> > ** **
>> >
>> > Regards,****
>> >
>> > ** **
>> >
>> > Hossein****
>> >
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucsd.edu/pipermail/ling-r-lang-l/attachments/20111010/543e01c5/attachment-0001.html
More information about the ling-r-lang-L
mailing list