Plotting mixed-effects model results with effects package

August 13, 2014
By

(This article was first published on Minding the Brain, and kindly contributed to R-bloggers)

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                                                     

> 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)
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 leave a comment for the author, please follow the link and comment on his blog: Minding the Brain.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.