Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
262 lines (174 sloc) 14.6 KB
---
title: Getting started with emmeans
author: Ariel Muldoon
date: '2019-03-25'
slug: getting-started-with-emmeans
categories:
- r
- statistics
tags:
- analysis
- teaching
- emmeans
draft: FALSE
description: "Post hoc comparisons are made easy in package emmeans. This post goes through some of the basics for those just getting started with the package."
---
```{r setup, echo = FALSE}
knitr::opts_chunk$set(comment = "#")
```
```{r, echo = FALSE}
# Make dataset here that will use later
# Two factors plus one log-normal response
set.seed(16)
dat = data.frame(f1 = rep(c("a", "c"), each = 10),
f2 = rep(c("1", "c"), times = 10) )
y = with(dat, 0 +
-1*(f1 == "a" & f2 == "c") +
-.5*(f1 == "c" & f2 == "1") +
1*(f1 == "c" & f2 == "c") +
rnorm(n = 20, mean = 0, sd = 1) )
dat$resp = round(exp(y), 1)
```
Package **emmeans** (formerly known as **lsmeans**) is enormously useful for folks wanting to do post hoc comparisons among groups after fitting a model. It has a very thorough set of vignettes (see the vignette topics [here](https://cran.r-project.org/web/packages/emmeans/vignettes/basics.html#more)), is very flexible with a ton of options, and works out of the box with a lot of different model objects (and can be extended to others `r emo::ji("+1")`).
I've started recommending **emmeans** all the time to students fitting models in R. However, I've found that often times students struggle a bit to get started using the package, possibly due to the sheer amount of flexibility and information in the vignettes.
I've put together some basic examples for using **emmeans**, meant to be a complement to the vignettes. Specifically this post will demonstrate a few of the built-in options for some standard post hoc comparisons; I will write a separate post about custom comparisons in **emmeans**.
*Disclaimer: This post is about using a package in R and so unfortunately does not focus on appropriate statistical practice for model fitting and post hoc comparisons.*
# R packages
I will load **magrittr** for the pipe in addition to **emmeans**.
```{r}
library(emmeans) # v. 1.3.3
library(magrittr) # v. 1.5
```
# The dataset and model
I've made a small dataset to use in this example. The response variable is `resp`, which comes from the log-normal distribution, and the two crossed factors of interest are `f1` and `f2`. Each factor has two levels: a control called `c` as well as a second non-control level.
```{r}
dat = structure(list(f1 = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("a",
"c"), class = "factor"), f2 = structure(c(1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("1",
"c"), class = "factor"), resp = c(1.6, 0.3, 3, 0.1, 3.2, 0.2,
0.4, 0.4, 2.8, 0.7, 3.8, 3, 0.3, 14.3, 1.2, 0.5, 1.1, 4.4, 0.4,
8.4)), row.names = c(NA, -20L), class = "data.frame")
str(dat)
```
The model I will use is a linear model with a log-transformed response and the two factors and their interaction as explanatory variables. This is the "true" model since I created these data so I'm skipping all model checks here (which I would not do in a real analysis).
Note I use `log(resp)` in the model rather than creating a new log-transformed variable. This will allow me to demonstrate one of the convenient options available in `emmeans()` later.
```{r}
fit1 = lm(log(resp) ~ f1 + f2 + f1:f2, data = dat)
```
# Built in comparisons with emmeans()
The **emmeans** package has helper functions for commonly used post hoc comparisons (aka *contrasts*). For example, we can do pairwise comparisons via `pairwise` or `revpairwise`, treatment vs control comparisons via `trt.vs.ctrl` or `trt.vs.ctrlk`, and even consecutive comparisons via `consec`.
The available built-in functions for doing comparisons are listed in the documentation for `?"contrast-methods"`.
# [All pairwise comparisons](#all-pairwise-comparisons)
One way to use `emmeans()`, which I use a lot, is to use formula coding for the comparisons. This formula is defined in the `specs` argument.
I will do all pairwise comparisons for all combinations of `f1` and `f2`. The built-in function `pairwise` is put on the left-hand side of the formula in `specs` and the factors with levels we want to compare among are on the right-hand side. Since I'm doing all pairwise comparisons, the combination of `f1` and `f2` are put in the formula.
The model object is passed to the first argument in `emmeans()`, `object`.
```{r}
emm1 = emmeans(fit1, specs = pairwise ~ f1:f2)
```
Using the formula in this way returns an object with two parts. The first part, called `emmeans`, is the estimated marginal means along with the standard errors and confidence intervals. We can pull these out with dollar sign notation, which I do below.
These results are all on the *model* scale, so in this case these are estimated mean log response for each `f1` and `f2` combination. Note the message that `emmeans()` gives us about results being on the log scale in the output. It knows the model is on the log scale because I used `log(resp)` as my response variable.
```{r}
emm1$emmeans
```
The second part of the output, called `contrasts`, contains the comparisons of interest. It is this section that we are generally most interested in when answering a question about differences among groups.
These results are also on the model scale (and we get the same message in this section), and [we'll want to put them on the original scale](#back-transforming-results).
The comparisons are accompanied by statistical tests of the null hypothesis of "no difference", but lack confidence interval (CI) limits by default. [We'll need to get these](#confidence-intervals-for-comparisons).
The `emmeans()` package automatically adjusts for multiple comparisons. Since we did all pairwise comparisons the package used a Tukey adjustment. [The type of adjustment can be changed](#changing-the-multiple-comparisons-adjustment).
```{r}
emm1$contrasts
```
# [Back-transforming results](#back-transforming-results)
Since I used a log transformation I can express the results as multiplicative differences in medians on the original (data) scale. We can always back-transform estimates and CI limits by hand, but in `emmeans()` we can use the `type` argument for this. Using `type = "response"` will return results on the original scale. This works when the transformation is explicit in the model (e.g., `log(resp)`) and works similarly for link functions in generalized linear models.
You'll see the message changes in the output once I do this, indicating things were back-transformed from the model scale. We also are reminded that the tests were done on the model scale.
In the `contrast` column in the `contrasts` section we can see the expression of the comparisons has changed from additive comparisons (via subtraction) shown above to multiplicative comparisons (via division).
```{r}
emmeans(fit1, specs = pairwise ~ f1:f2, type = "response")
```
# [Changing the multiple comparisons adjustment](#changing-the-multiple-comparisons-adjustment)
The `adjust` argument can be used to change the type of multiple comparisons adjustment. All available options are listed and described in the documentation for `summary.emmGrid` under the section *P-value adjustments*.
One option is to skip multiple comparisons adjustments all together, using `adjust = "none"`. If we use this the message about multiple comparisons disappears (since we didn't use one).
```{r}
emm1.1 = emmeans(fit1, specs = pairwise ~ f1:f2, type = "response", adjust = "none")
emm1.1
```
# [Confidence intervals for comparisons](#confidence-intervals-for-comparisons)
We will almost invariably want to report confidence intervals for any comparisons of interest. We need a separate function to get these. Here is an example using the `confint()` function with the default 95% CI (the confidence level can be changed, see `?confint.emmGrid`). I use the pipe to pass the `contrasts` into the `confint()` function.
```{r}
emm1.1$contrasts %>%
confint()
```
The `confint()` function returns confidence intervals but gets rid of the statistical tests. Some people will want to also report the test statistics and p-values. In this case, we can use `summary()` instead of `confint()`, with `infer = TRUE`.
```{r}
emm1.1$contrasts %>%
summary(infer = TRUE)
```
# [Putting results in a data.frame](#putting-results-in-a-data.frame)
One of the really nice things about `emmeans()` is that it makes it easy to get the results into a nice format for making tables or graphics of results. This is because the results can be converted into a data.frame via `as.data.frame()`.
This is true for both `confint()` and `summary()`, although I'll just show `confint()` here. After getting the confidence intervals I pipe the results into `as.data.frame()`. I assign this object a name so I can use it for, e.g., making a graph of results.
```{r}
emm1.1_contrasts = emm1.1$contrasts %>%
confint() %>%
as.data.frame()
emm1.1_contrasts
```
If needed, the estimated marginal means can also be put into a data.frame.
```{r}
emm1.1$emmeans %>%
as.data.frame()
```
# [Within group comparisons](#within-group-comparisons)
While we *can* do all pairwise comparisons, there are certainly plenty of situations where the research question dictates that we only want a specific set of comparisons. A common example of this is when we want to compare the levels of one factor within the levels of another. Here I'll show comparisons among levels of `f1` for each level of `f2`.
The only thing that changes is the right-hand side of the `specs` formula. The code `f1|f2` translates to "compare levels of `f1` within each level of `f2`".
```{r}
emm2 = emmeans(fit1, specs = pairwise ~ f1|f2, type = "response")
emm2
```
You can see there is no message about a multiple comparisons adjustment in the above set of comparisons. This is because the package default is to correct for the number of comparisons *within* each group instead of across groups. In this case there is only a single comparison in each group.
If we consider the family of comparisons to be all comparisons regardless of group and want to correct for multiple comparisons, we can do so via `rbind.emmGrid`.
Here is an example of passing `contrasts` to `rbind()` to correct for multiple comparisons. The default adjustment is Bonferroni, which can be much too conservative when the number of comparisons is large. You can control the multiple comparisons procedure via `adjust`.
The results of `rbind()` can also be used with `summary()`, `confint()`, and/or `as.data.frame()`.
```{r}
emm2$contrasts %>%
rbind()
```
# [Main effects comparisons](#main-effects-comparisons)
Even if we have multiple factors in the model, complete with an interaction term, we can still do "overall" comparisons among groups if our research question indicated that main effects were an important thing to estimate.
Doing main effects in the presence of an interaction means we *average over* the levels of the other factor(s). The `emmeans()` function gives both a warning about the interaction and a message indicating which factor was averaged over to remind us of this.
Here is the estimated main effect of `f1`. Since we are only interested in overall comparisons of that factor it is the only factor given on the right-hand side of the `specs` formula.
```{r}
emmeans(fit1, specs = pairwise ~ f1)
```
# [Treatment vs control example](#treatment-vs-control-example)
The **emmeans** package has built-in helper functions for comparing each group mean to the control mean. If the control group is the in the first row of the `emmeans` section of the output, this set of comparisons can be requested via `trt.vs.ctrl`.
Note the default multiple comparisons adjustment is a Dunnett adjustment.
```{r}
emmeans(fit1, specs = trt.vs.ctrl ~ f1:f2)
```
Using `trt.vs.ctrl` means we ended up comparing each group mean to the "a 1" group since it is in the first row. In the example I'm using the control group, "c c", is actually the *last* group listed in the `emmeans` section. When the control group is the last group in `emmeans` we can use `trt.vs.ctrlk` to get the correct set of comparisons.
```{r}
emmeans(fit1, specs = trt.vs.ctrlk ~ f1:f2)
```
That gives us what we want in this case. However, if the control group was some other group, like "c 1", we could use `trt.vs.ctrlk` with the `ref` argument to define which row in the `emmeans` section represents the control group.
The "c 1" group is the second row in the `emmeans` so we can use `ref = 2` to define this group as the control group.
```{r}
emmeans(fit1, specs = trt.vs.ctrlk ~ f1:f2, ref = 2)
```
Finally, if we want to reverse the order of subtraction in the treatment vs control comparisons we can use the `reverse` argument.
```{r}
emmeans(fit1, specs = trt.vs.ctrlk ~ f1:f2, ref = 2, reverse = TRUE)
```
# [Alternative code for comparisons](#alternative-code-for-comparisons)
The `emmeans()` package offers the option to do comparisons in two steps instead of in one step the way I have been using it as above. I personally find this alternative most useful when doing custom comparisons, and I think it's good to introduce it now so it looks familiar. This alternative keeps the estimated marginal means and the comparisons of interest in separate objects, which may be attractive in some situations.
The first step is to use `emmeans()` to calculate the marginal means of interest. We still use the formula in `specs` with the factor(s) of interest on the right-hand side but no longer put anything on the left-hand side of the tilde.
We can still use `type` in `emmeans()` but cannot use `adjust` (since we don't adjust for multiple comparisons until we've actually done comparisons `r emo::ji("wink")`).
```{r}
emm3 = emmeans(fit1, specs = ~ f1:f2, type = "response")
emm3
```
We then get the comparisons we want in a second step using the `contrast()` function. We request the comparisons we want via `method`. When using built-in comparisons like I am here, we give the comparison function name as a string (meaning in quotes). Also see the `pairs()` function, which is for the special case of all pairwise comparisons.
We can use `adjust` in `contrast()` to change the multiple comparisons adjustment.
```{r}
contrast(emm3, method = "pairwise", adjust = "none")
```
We can follow the `contrast()` argument with `summary()` or `confint()` and then put the results into a data.frame for plotting/saving. Again, I think the real strength of `contrast()` comes when we want custom comparisons, and I'll demonstrate these in my next post.
You can’t perform that action at this time.