Plotting mixed-effects model results with effects package
[This article was first published on Minding the Brain, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
As separate by-subjects and by-items analyses have been replaced by mixed-effects models with crossed random effects of subjects and items, I’ve often found myself wondering about the best way to plot data. The simple-minded means and SE from trial-level data will be inaccurate because they won’t take the nesting into account. If I compute subject means and plot those with by-subject SE, then I’m plotting something different from what I analyzed, which is not always terrible, but definitely not ideal. It seems intuitive that the condition means and SE’s are computable from the model’s parameter estimates, but that computation is not trivial, particularly when you’re dealing with interactions. Or, rather, that computation was not trivial until I discovered the effects package.
To show how this would work, I pillaged some data from a word-to-picture matching pilot study. Younger adults (college students) and older adults (mostly 50s and 60s) did a word-to-picture matching task in the presence of either cohort competitors (camera – camel) or semantic competitors (lion – tiger).
> summary(RT.demo)
Subject Target Condition ACC RT Group
2102 : 40 chicken: 36 Cohort :690 Min. :1 Min. :1503 YC:701
2103 : 40 hat : 36 Semantic:675 1st Qu.:1 1st Qu.:2131 OC:664
2104 : 40 penny : 36 Median :1 Median :2362
2106 : 40 potato : 36 Mean :1 Mean :2442
2109 : 40 radio : 36 3rd Qu.:1 3rd Qu.:2684
2116 : 40 stool : 36 Max. :1 Max. :4847
(Other):1125 (Other):1149
Or for people who like dynamite plots:
> ggplot(x, aes(Group, fit, color=Condition, fill=Condition)) + geom_bar(stat=”identity”, position=”dodge”) + geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4, position=position_dodge(width=0.9)) + theme_bw(base_size=12)
To show how this would work, I pillaged some data from a word-to-picture matching pilot study. Younger adults (college students) and older adults (mostly 50s and 60s) did a word-to-picture matching task in the presence of either cohort competitors (camera – camel) or semantic competitors (lion – tiger).
> summary(RT.demo)
Subject Target Condition ACC RT Group
2102 : 40 chicken: 36 Cohort :690 Min. :1 Min. :1503 YC:701
2103 : 40 hat : 36 Semantic:675 1st Qu.:1 1st Qu.:2131 OC:664
2104 : 40 penny : 36 Median :1 Median :2362
2106 : 40 potato : 36 Mean :1 Mean :2442
2109 : 40 radio : 36 3rd Qu.:1 3rd Qu.:2684
2116 : 40 stool : 36 Max. :1 Max. :4847
(Other):1125 (Other):1149
> ggplot(RT.demo, aes(Condition, RT, fill=Group, color=Group)) + geom_violin() + theme_bw(base_size=12)
Not surprisingly, the response times for older adults are slower than for younger adults, but it looks like this might be particularly true in the presence of semantic competitors. Let’s test that with a mixed model with crossed random effects of subjects and items.
> m <- lmer(RT ~ Condition*Group + (Condition | Subject) + (1 | Target), data=RT.demo)
> coef(summary(m))
Estimate Std. Error t value
(Intercept) 2230.057 64.749 34.44
ConditionSemantic -7.881 68.565 -0.11
GroupOC 413.287 65.097 6.35
ConditionSemantic:GroupOC 104.096 33.110 3.14
So it looks like the older adults are about 400ms slower than the younger adults in the cohort condition and another 100ms slower in the semantic condition. Now we can use the effects package to convert these parameter estimates into condition mean and SE estimates. The key function is effect(), which takes a term from the model and the model object. We can use summary() on the effect list object to get the information we need.
> library(effects)
> ef <- effect("Condition:Group", m)
> summary(ef)
Condition*Group effect
Group
Condition YC OC
Cohort 2230.057 2643.344
Semantic 2222.176 2739.559
Lower 95 Percent Confidence Limits
Group
Condition YC OC
Cohort 2103.037 2516.026
Semantic 2088.161 2605.384
Upper 95 Percent Confidence Limits
Group
Condition YC OC
Cohort 2357.076 2770.662
Semantic 2356.190 2873.734
For the purposes of plotting, we want to convert the effect list object into a data frame. Conveniently, there is a as.data.frame() function:
> x <- as.data.frame(ef)
> x
Condition Group fit se lower upper
1 Cohort YC 2230.057 64.74945 2103.037 2357.076
2 Semantic YC 2222.176 68.31514 2088.161 2356.190
3 Cohort OC 2643.344 64.90167 2516.026 2770.662
4 Semantic OC 2739.559 68.39711 2605.384 2873.734
Now we can plot this:
> ggplot(x, aes(Condition, fit, color=Group)) + geom_point() + geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4) + theme_bw(base_size=12)
> ggplot(x, aes(Group, fit, color=Condition, fill=Condition)) + geom_bar(stat=”identity”, position=”dodge”) + geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4, position=position_dodge(width=0.9)) + theme_bw(base_size=12)
To leave a comment for the author, please follow the link and comment on their blog: Minding the Brain.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.