[R-lang] Simple effects in mixed logit models
T. Florian Jaeger
tiflo at csli.stanford.edu
Fri Dec 18 16:44:12 PST 2009
Hi Roger,
I have a quite lot of missing responses, and as you can see in the summary
> of the results, there are correlations between the “main” effects and the
> interaction. But you say “not the fixed effect correlation”, so are you
> referring to a different correlation?
>
I just meant that the specific numbers I gave were those of the correlations
between the two variables (cor()). Fixed effect correlations are a somewhat
complicated beast, so I wasn't talking about them, but usually when the
correlations are high, the fixed effect correlations are high, too.
> I haven’t done the model comparisons yet. You said:
>
> In any case, since ANOVA are omnibus test, you can get at least that by
> simple comparing the above model to a model without the interaction or
> without any of the simple (main) effects (using anova(model1, model2)).
> Model comparison is robust against collinearity.
>
>
>
> I find it hard to conceptualize what the model comparisons tell me. What
> can I conclude if the model with which I got the simple effect differs from
> a model without the interaction or without any of the simple effects?
>
I had a somewhat hard time finding what you did in your code below, but I
was thinking that you can just compare the models with and without a simple
effect. if the comparison is significant, the simple effect matters. I put
up some code that might be generally helpful at
http://hlplab.wordpress.com/2009/12/18/some-r-code-to-understand-the-difference-between-treatment-and-sum-anova-style-coding/,
but I didn't have time to clean things up. Let me know what you think.
>
>
> My previous analyses were done with R 2.2.0. I have now updated to 2.10.1
> and suddenly, all the coefficients are different, which I find a bit
> worrying. The p-values are now quite a bit larger, although the new R
> version doesn’t change any of the conclusions that I can draw from the data.
>
wow. That is worrisome. I'd be interested to hear what other people's
experience with R 2.10.1 is. I haven't upgraded yet.
Btw, I had a quick look at your code below and I have two thoughts:
- you can treatment code automatically, by saying:
- contrasts(factor) <- contr.treatment(numer_of_levels_of_factor)
- e.g. contrasts(language) <- contr.treatment(2)
- the contrasts() command creates attributes to the variable which
regression functions are sensitive to. i.e. if you put a factor into a
regression function (such as lmer) the regression will use the
contrast-coding specified by contrasts. to see what the contrasts are,
simply enter contrasts(factor).
- as a matter of fact, treatment coding is the default in R. so, if
you just throw in a factor into a regression, R will treatment code it.
- an easy way to create a new variable that contains the contrast coding
uses the ifelse statement:
- data$LanguageSwedish <- ifelse(data$language == "swe", 1, 0)
maybe useful?
Florian
>
>
> Roger
>
>
>
>
>
> > library(lme4)
>
> Loading required package: Matrix
>
> Loading required package: lattice
>
> > L2prim = read.table("forLME.txt", header=TRUE) # open data file
>
> > head(L2prim) # show first 6 lines
>
> LIST PID ITEM language prime score
>
> 1 1 2 4 swe do 1
>
> 2 1 2 7 swe do 0
>
> 3 1 2 8 swe pp 0
>
> 4 1 2 20 swe pp 0
>
> 5 1 2 22 eng pp 0
>
> 6 1 2 23 swe pp 0
>
> >
>
> > # effect code vars so that they can be interpreted as main effects in
> ANOVAs
>
> > L2prim$langcent <- scale(as.numeric(L2prim$language)) # effect code
> language
>
> > L2prim$primecent <- scale(as.numeric(L2prim$prime)) # effect code prime
>
> >
>
> > L2prim$PID=as.factor(L2prim$PID) # convert subject no to factor
>
> > L2prim$ITEM=as.factor(L2prim$ITEM) # convert item no to factor
>
> > head(L2prim)
>
> LIST PID ITEM language prime score langcent primecent
>
> 1 1 2 4 swe do 1 1.008740 -0.9742244
>
> 2 1 2 7 swe do 0 1.008740 -0.9742244
>
> 3 1 2 8 swe pp 0 1.008740 1.0248334
>
> 4 1 2 20 swe pp 0 1.008740 1.0248334
>
> 5 1 2 22 eng pp 0 -0.989767 1.0248334
>
> 6 1 2 23 swe pp 0 1.008740 1.0248334
>
> >
>
> > # mixed logit model with variables effect coded (main effects as in
> ANOVA)
>
> > primanal=lmer(score ~langcent*primecent + (1|PID) + (1|ITEM),
> data=L2prim, family = "binomial")
>
> > summary(primanal)
>
> Generalized linear mixed model fit by the Laplace approximation
>
> Formula: score ~ langcent * primecent + (1 | PID) + (1 | ITEM)
>
> Data: L2prim
>
> AIC BIC logLik deviance
>
> 650.5 677.2 -319.2 638.5
>
> Random effects:
>
> Groups Name Variance Std.Dev.
>
> ITEM (Intercept) 1.9705 1.4037
>
> PID (Intercept) 2.2791 1.5097
>
> Number of obs: 632, groups: ITEM, 40; PID, 32
>
>
>
> Fixed effects:
>
> Estimate Std. Error z value Pr(>|z|)
>
> (Intercept) -1.06339 0.37331 -2.849 0.004392 **
>
> langcent 0.01506 0.10822 0.139 0.889334
>
> primecent 0.36986 0.11024 3.355 0.000793 ***
>
> langcent:primecent -0.10630 0.10898 -0.975 0.329365
>
> ---
>
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
>
> Correlation of Fixed Effects:
>
> (Intr) lngcnt prmcnt
>
> langcent -0.001
>
> primecent -0.028 -0.030
>
> lngcnt:prmc 0.008 -0.055 -0.008
>
> >
>
> >
>
> > # Analysis of simple effect of prime for Swedish (coded as 0)
>
> >
>
> > L2prim[,2:3] # Outputs col 2 and 3 (to check no of rows)
>
> PID ITEM
>
> 1 2 4
>
> 2 2 7
>
> 3 2 8
>
> 4 2 20
>
> 5 2 22
>
> 6 2 23
>
>
>
> 629 34 35
>
> 630 34 36
>
> 631 34 37
>
> 632 34 38
>
> >
>
> > # create column for language that is treatment coded
>
> > x =c(1:632) # creates 632 rows of numbers
>
> > x[1:632] = 0 # sets the numbers to 0
>
> > cbind(L2prim,x) # adds the 0s to L2prim
>
> LIST PID ITEM language prime score langcent primecent x
>
> 1 1 2 4 swe do 1 1.008740 -0.9742244 0
>
> 2 1 2 7 swe do 0 1.008740 -0.9742244 0
>
> 3 1 2 8 swe pp 0 1.008740 1.0248334 0
>
> 4 1 2 20 swe pp 0 1.008740 1.0248334 0
>
> 5 1 2 22 eng pp 0 -0.989767 1.0248334 0
>
> 6 1 2 23 swe pp 0 1.008740 1.0248334 0
>
>
>
> 629 4 34 35 eng pp 1 -0.989767 1.0248334 0
>
> 630 4 34 36 swe do 1 1.008740 -0.9742244 0
>
> 631 4 34 37 eng pp 1 -0.989767 1.0248334 0
>
> 632 4 34 38 eng do 0 -0.989767 -0.9742244 0
>
> > L2prim[L2prim[,4] =="swe",9] = 0 # recodes swe (col 4) into 1 (col 9)
>
> > L2prim[L2prim[,4] =="eng",9] = 1 # recodes eng (col 4) into 0 (col 9)
>
> >
>
> > # create column for prime that is treatment coded
>
> > y =c(1:632) # creates 632 rows of numbers
>
> > y[1:632] = 0 # sets the numbers to 0
>
> > cbind(L2prim,y) # adds the 0s to L2prim
>
> LIST PID ITEM language prime score langcent primecent V9 y
>
> 1 1 2 4 swe do 1 1.008740 -0.9742244 0 0
>
> 2 1 2 7 swe do 0 1.008740 -0.9742244 0 0
>
> 3 1 2 8 swe pp 0 1.008740 1.0248334 0 0
>
> 4 1 2 20 swe pp 0 1.008740 1.0248334 0 0
>
> 5 1 2 22 eng pp 0 -0.989767 1.0248334 1 0
>
> 6 1 2 23 swe pp 0 1.008740 1.0248334 0 0
>
>
>
> 629 4 34 35 eng pp 1 -0.989767 1.0248334 1 0
>
> 630 4 34 36 swe do 1 1.008740 -0.9742244 0 0
>
> 631 4 34 37 eng pp 1 -0.989767 1.0248334 1 0
>
> 632 4 34 38 eng do 0 -0.989767 -0.9742244 1 0
>
> > L2prim[L2prim[,5] =="do",10] = 1 # recodes do (col 5) into 1 (col 10)
>
> > L2prim[L2prim[,5] =="pp",10] = 0 # recodes pp (col 5) into 0 (col 10)
>
> >
>
> > L2prim$langtreat=as.factor(L2prim$V9) # convert language (v9) into factor
>
> > L2prim$primetreat=as.factor(L2prim$V10) # convert prime (v10) into factor
>
> > L2prim[1:20,] # Outputs rows 1-20
>
> LIST PID ITEM language prime score langcent primecent V9 V10
> langtreat
>
> 1 1 2 4 swe do 1 1.008740 -0.9742244 0 1
> 0
>
> 2 1 2 7 swe do 0 1.008740 -0.9742244 0 1
> 0
>
> 3 1 2 8 swe pp 0 1.008740 1.0248334 0 0
> 0
>
> 4 1 2 20 swe pp 0 1.008740 1.0248334 0 0
> 0
>
> 5 1 2 22 eng pp 0 -0.989767 1.0248334 1 0
> 1
>
> 6 1 2 23 swe pp 0 1.008740 1.0248334 0 0
> 0
>
>
>
> 20 1 3 10 eng pp 0 -0.989767 1.0248334 1 0
> 1
>
> primetreat
>
> 1 1
>
> 2 1
>
> 3 0
>
> 4 0
>
> 5 0
>
> 6 0
>
>
>
> 20 0
>
> >
>
> > # mixed logit model with variables treatment coded
>
> > primanal=lmer(score ~langtreat*primetreat + (1|PID) + (1|ITEM),
> data=L2prim, family = "binomial")
>
> > summary(primanal)
>
> Generalized linear mixed model fit by the Laplace approximation
>
> Formula: score ~ langtreat * primetreat + (1 | PID) + (1 | ITEM)
>
> Data: L2prim
>
> AIC BIC logLik deviance
>
> 650.5 677.2 -319.2 638.5
>
> Random effects:
>
> Groups Name Variance Std.Dev.
>
> ITEM (Intercept) 1.9705 1.4038
>
> PID (Intercept) 2.2791 1.5097
>
> Number of obs: 632, groups: ITEM, 40; PID, 32
>
>
>
> Fixed effects:
>
> Estimate Std. Error z value Pr(>|z|)
>
> (Intercept) -0.7790 0.4155 -1.875 0.0608 .
>
> langtreat1 0.1875 0.3022 0.621 0.5348
>
> primetreat1 -0.5251 0.3100 -1.694 0.0903 .
>
> langtreat1:primetreat1 -0.4246 0.4354 -0.975 0.3294
>
> ---
>
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
>
> Correlation of Fixed Effects:
>
> (Intr) lngtr1 prmtr1
>
> langtreat1 -0.364
>
> primetreat1 -0.354 0.477
>
> lngtrt1:pr1 0.262 -0.699 -0.703
>
> > # primetreat is the simple effect of prime on Swedish
>
> > # the sum of the estimates for primetreat and the interaction (0.9497) is
> the simple effect
>
> > # of prime on English (no z and p values can be calculated)
>
> >
>
> >
>
> > # Analysis of simple effect of prime for English (coded as 0)
>
> >
>
> > L2prim[,2:3] # Outputs col 2 and 3 (to check no of rows)
>
> PID ITEM
>
> 1 2 4
>
> 2 2 7
>
> 3 2 8
>
> 4 2 20
>
> 5 2 22
>
> 6 2 23
>
>
>
> 629 34 35
>
> 630 34 36
>
> 631 34 37
>
> 632 34 38
>
> >
>
> > x =c(1:632) # creates 632 rows of numbers
>
> > x[1:632] = 0 # sets the numbers to 0
>
> > cbind(L2prim,x) # adds the 0s to L2prim
>
> LIST PID ITEM language prime score langcent primecent V9 V10
> langtreat
>
> 1 1 2 4 swe do 1 1.008740 -0.9742244 0 1
> 0
>
> 2 1 2 7 swe do 0 1.008740 -0.9742244 0 1
> 0
>
> 3 1 2 8 swe pp 0 1.008740 1.0248334 0 0
> 0
>
> 4 1 2 20 swe pp 0 1.008740 1.0248334 0 0
> 0
>
> 5 1 2 22 eng pp 0 -0.989767 1.0248334 1 0
> 1
>
> 6 1 2 23 swe pp 0 1.008740 1.0248334 0 0
> 0
>
>
>
> 629 4 34 35 eng pp 1 -0.989767 1.0248334 1 0
> 1
>
> 630 4 34 36 swe do 1 1.008740 -0.9742244 0 1
> 0
>
> 631 4 34 37 eng pp 1 -0.989767 1.0248334 1 0
> 1
>
> 632 4 34 38 eng do 0 -0.989767 -0.9742244 1 1
> 1
>
> primetreat x
>
> 1 1 0
>
> 2 1 0
>
> 3 0 0
>
> 4 0 0
>
> 5 0 0
>
> 6 0 0
>
>
>
> 629 0 0
>
> 630 1 0
>
> 631 0 0
>
> 632 1 0
>
> > L2prim[L2prim[,4] =="swe",9] = 1 # recodes swe (col 4) into 1 (col 9)
>
> > L2prim[L2prim[,4] =="eng",9] = 0 # recodes eng (col 4) into 0 (col 9)
>
> >
>
> > y =c(1:632) # creates 632 rows of numbers
>
> > y[1:632] = 0 # sets the numbers to 0
>
> > cbind(L2prim,y) # adds the 0s to L2prim
>
> LIST PID ITEM language prime score langcent primecent V9 V10
> langtreat
>
> 1 1 2 4 swe do 1 1.008740 -0.9742244 1 1
> 0
>
> 2 1 2 7 swe do 0 1.008740 -0.9742244 1 1
> 0
>
> 3 1 2 8 swe pp 0 1.008740 1.0248334 1 0
> 0
>
> 4 1 2 20 swe pp 0 1.008740 1.0248334 1 0
> 0
>
> 5 1 2 22 eng pp 0 -0.989767 1.0248334 0 0
> 1
>
> 6 1 2 23 swe pp 0 1.008740 1.0248334 1 0
> 0
>
>
>
> 629 4 34 35 eng pp 1 -0.989767 1.0248334 0 0
> 1
>
> 630 4 34 36 swe do 1 1.008740 -0.9742244 1 1
> 0
>
> 631 4 34 37 eng pp 1 -0.989767 1.0248334 0 0
> 1
>
> 632 4 34 38 eng do 0 -0.989767 -0.9742244 0 1
> 1
>
> primetreat y
>
> 1 1 0
>
> 2 1 0
>
> 3 0 0
>
> 4 0 0
>
> 5 0 0
>
> 6 0 0
>
>
>
> 628 0 0
>
> 629 0 0
>
> 630 1 0
>
> 631 0 0
>
> 632 1 0
>
> > L2prim[L2prim[,5] =="do",10] = 1 # recodes do (col 5) into 1 (col 10)
>
> > L2prim[L2prim[,5] =="pp",10] = 0 # recodes pp (col 5) into 0 (col 10)
>
> >
>
> > L2prim$langtreat=as.factor(L2prim$V9)
>
> > L2prim$primetreat=as.factor(L2prim$V10)
>
> > L2prim[1:20,] # Outputs rows 1-20
>
> LIST PID ITEM language prime score langcent primecent V9 V10
> langtreat
>
> 1 1 2 4 swe do 1 1.008740 -0.9742244 1 1
> 1
>
> 2 1 2 7 swe do 0 1.008740 -0.9742244 1 1
> 1
>
> 3 1 2 8 swe pp 0 1.008740 1.0248334 1 0
> 1
>
> 4 1 2 20 swe pp 0 1.008740 1.0248334 1 0
> 1
>
> 5 1 2 22 eng pp 0 -0.989767 1.0248334 0 0
> 0
>
> 6 1 2 23 swe pp 0 1.008740 1.0248334 1 0
> 1
>
>
>
> 20 1 3 10 eng pp 0 -0.989767 1.0248334 0 0
> 0
>
> primetreat
>
> 1 1
>
> 2 1
>
> 3 0
>
> 4 0
>
> 5 0
>
> 6 0
>
>
>
> 20 0
>
> >
>
> > primanal=lmer(score ~langtreat*primetreat + (1|PID) + (1|ITEM),
> data=L2prim, family = "binomial")
>
> > summary(primanal)
>
> Generalized linear mixed model fit by the Laplace approximation
>
> Formula: score ~ langtreat * primetreat + (1 | PID) + (1 | ITEM)
>
> Data: L2prim
>
> AIC BIC logLik deviance
>
> 650.5 677.2 -319.2 638.5
>
> Random effects:
>
> Groups Name Variance Std.Dev.
>
> ITEM (Intercept) 1.9705 1.4038
>
> PID (Intercept) 2.2791 1.5097
>
> Number of obs: 632, groups: ITEM, 40; PID, 32
>
>
>
> Fixed effects:
>
> Estimate Std. Error z value Pr(>|z|)
>
> (Intercept) -0.5915 0.4153 -1.424 0.15435
>
> langtreat1 -0.1875 0.3022 -0.621 0.53489
>
> primetreat1 -0.9496 0.3095 -3.068 0.00215 **
>
> langtreat1:primetreat1 0.4244 0.4354 0.975 0.32973
>
> ---
>
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
>
> Correlation of Fixed Effects:
>
> (Intr) lngtr1 prmtr1
>
> langtreat1 -0.363
>
> primetreat1 -0.354 0.507
>
> lngtrt1:pr1 0.247 -0.699 -0.702
>
> > # primetreat is the simple effect of prime on English
>
> > # the sum of the estimates for primetreat and the interaction (0.5252) is
> the simple effect
>
> > # of prime on Swedish (no z and p values can be calculated)
>
> The University of Dundee is a registered Scottish charity, No: SC015096
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://pidgin.ucsd.edu/pipermail/r-lang/attachments/20091219/d8d821dd/attachment.htm>
More information about the R-lang
mailing list