Overview

Here we will be giving an overview of estimated marginal means and how to calculate them in R. Most of the information uses examples from Getting started with emmeans. Additionally, review on the emmeans() function in R is available at rvlenth/emmeans

Background

According to ChatGPT, estimated marginal means represent an average response (mean) of a dependent variable at a certain level of a predictor while accounting for the other covariates in the model, or holding them constant. Thus, estimate marginal means are the expected value of Y in models that have more than one predictor:

$\hat{Y}_{marginal} = \beta_{0} + \beta_{1}X_{1} + \beta_{2}X_{2} + ... + \beta_{k}X_{k}$

Where:

• $\hat{Y}_{marginal}$ is the estimated marginal mean for a certain level of the factor or covariate of interest
• $\beta_{0}$ is the intercept, and $\beta_{1}, \beta_{2},..., \beta_{k}$ are the coefficients of the model.
• $X_{1}, X_{2},...,X_{k}$ are the values of the other covariates that are held in the calculation.

Essentially, emmeans() provides the mean outcome for each factor level, after adjusting for the other variables in the model.

Additionally, there is another function named emtrends to be familiar with. While emmeans() is calculated the estimated group mean of levels in a factor, emtrends is calculating estimated slopes or rate of change! In other words, it is adjusting slopes for continuous predictors at each level of a factor.

Getting Started

Creating dummy data with one outcome (resp) and two factors (f1 and f2).

# Creating dummy data
dat = data.frame(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),
                 f1 = factor(c("a","a","a","a","a",
                               "a","a","a","a","a","c","c","c","c","c",
                               "c","c","c","c","c")),
                 f2 = factor(c("1","c","1","c","1",
                               "c","1","c","1","c","1","c","1","c","1",
                               "c","1","c","1","c")))

Creating a linear model with this data while doing a log transformation on the outcome.

fit1 = lm(log(resp) ~ f1 + f2 + f1:f2, data = dat)
summary(fit1)

Model Output

• By looking at this output, we can see that there is no main effect of f1, a main effect of f2, and a significant f1 and f2 interaction.

Different comparisons to be made

Since emmeans() is comparing the means of the levels of a factor, there are ways to modify how those means should be calculated or how the comparisons between the means should be done. There is a specs argument within emmeans() where any of the following parameters can be used:

• pairwise
• revpairwise
• trt.vs.vtrl
• trt.vs.ctrlk
• consec

All pairwise comparisons

Using pairwise as a parameter for the interaction will produce a comparison between all possible level combinations in a pair manner.

emm1 = emmeans(fit1, specs = pairwise ~ f1:f2)

emmeans() output

• When we run emmeans() like this, the output produces a list with two datasets. The first named 'emmeans' contains the estimated marginal means of each combination of levels from the factors. In other words, we can get the expected mean for a specific group defined by one level of one factor and one level from another factor. These estimated means are in the same scale unit as the outcome was when ran in the model! So in this example, these estimated marginal means are all the log of Y.

• The other dataset is named 'contrast', which is the more important part. Here we have differences between different combinations of the estimated marginal means calculated in the dataset above. Along with the difference, we also have p-values associated with them, indicating if those two estimated marginal means are significantly different from each other!

Back-transforming results

• Since we decided to log transform our outcome in the model, we can revert this back so that the produced output from emmeans() is in the regular/original scale units of the outcome. To do this we just add the type = "response".

emmeans(fit1, specs = pairwise ~ f1:f2, type = "response")

emmeans() output

• As mentioned, all of the math is done in the log-transformed scale. Here those estimated marginal means are converted back into the original scale by exponentiating them (log-inverse). Since the contrasts were being subtracted in the log-transformed scale, now in the original they are ratios. These ratio represent the mean of group A divided by the mean of group B. Thus, if there is a ratio of 1.5, then that means the mean of Group A is 50% higher than the mean of group B. If the ratio is .8 for example, then that means the mean of group A is (1-.8 = .2) 20% lower than the mean of Group B.

Changing the multiple comparison adjustment

The emmeans() function allows you to choose how to adjust the p-values to account for the increased likelihood of Type I errors. You can just introduce the 'adjust' argument and make it equal to one of the following:

• "none"
• "tukey"
• "bonferroni"
• "sidak"
• "holm"
• "fdr"
• "scheffe"

You can learn more by writing this into R

?summary.emmGrid

Confidence Intervals

• If you want confidence intervals with your output- you can take the emmeans object and use the confint() function on it or a better approach is to use the summary() function on it while having the parameter 'infer = TRUE'. This will provide both the confidence interval and the p-value associated with the contrast.

summary(emm2, infer = TRUE)

Within group comparisons

While we can do pairwise comparisons for all levels of each predictor- we could also just do comparisons of the levels of one factor within the levels of another. The code below is saying compare the level of 'f1' for each level of 'f2'. Since 'f1' only has two levels, no correction for multiple comparison is implemented- this will also be obvious since there will be no message on this too.

emm2 = emmeans(fit1, specs = pairwise ~ f1|f2, type = "response")
emm2

BE CAREFUL, if you use the rbind() function to tidy your contrasts- you will automatically trigger a bonferroni correction to your p-values. If you do not want this to occur, uses data.frame() instead. Why does this happen? Emmeans assumes you want to compare differences not just stratified within levels but now across levels- so this requires another correction! (According to ChatGPT).

# BAD (I think)
rbind(emm2$contrasts)

# Good
data.frame(emm2$contrasts)

Main effects comparisons

Even if there is a significant interaction term, we can still investigate main effects if interested using emmeans(). However, the output will give you a warning to interpret the effects with caution due to interactions effects in the model.

emmeans(fit1, specs = pairwise ~ f1)

Treatment vs Control Contrasts

There are default contrast settings that you can use so your pairwise comparisons are more meaningful. For example, if the first group is a control group, and you only need comparisons to be against the control group- then we can specify this using 'specs = trt.vs.ctrl'.

emmeans(fit1, specs = trt.vs.ctrl ~ f1:f2)

emmeans() output

• Notice all of the estimated means are displayed at the top and the contrasts in the bottom. The contrasts are all comparisons against 'a 1', which is the first group in the $emmeans part. If the control group is not in the first row, then we can modify our code to change the comparison to a different row by using 'specs = trt.vs.ctrlk' and 'ref = 2', with the latter specifying which row the control group is in.

emmeans(fit1, specs = trt.vs.ctrlk ~ f1:f2, ref = 2)

• We can also reverse the order of subtraction by using the reverse argument.

emmeans(fit1, specs = trt.vs.ctrlk ~ f1:f2, ref = 2, reverse = TRUE)

Alternative Code for comparison

• So far, we have been using the arguments within the emmeans() function to create the multiple comparisons ('specs = pairwise' or 'specs = trt.vs.ctrlk'). This is perfectly fine for simple mean pairwise comparisons. However, we can also do these comparisons outside of emmeans() by using either contrast() or pairs(). Thus, we can create an object using emmeans but leave out the 'specs' portion entirely, only generating the estimated marginal means and not the contrasts.

emm3 = emmeans(fit1, specs = ~ f1:f2, type = "response")
emm3

Next we can use contrast() and specify the method of comparison as "pairwise". We can also pick what type of p-value adjustment we want if any.

contrast(emm3, method = "pairwise", adjust = "none")

• Why is this important? We use this approach when we want to do custom comparisons! View the next wiki page for more information.