Basics of estimated marginal means

The point

The point is that the marginal means of cell.means give equal weight to each cell. In many situations (especially with experimental data), that is a much fairer way to compute marginal means, in that they are not biased by imbalances in the data. We are, in a sense, estimating what the marginal means would be, had the experiment been balanced. Estimated marginal means (EMMs) serve that need.

All this said, there are certainly situations where equal weighting is not appropriate. Suppose, for example, we have data on sales of a product given different packaging and features. The data could be unbalanced because customers are more attracted to some combinations than others. If our goal is to understand scientifically what packaging and features are inherently more profitable, then equally weighted EMMs may be appropriate; but if our goal is to predict or maximize profit, the ordinary marginal means provide better estimates of what we can expect in the marketplace.

Back to Contents

Model and reference grid

Estimated marginal means are based on a model – not directly on data. The basis for them is what we call the reference grid for a given model. To obtain the reference grid, consider all the predictors in the model. Here are the default rules for constructing the reference grid

• For each predictor that is a factor, use its levels (dropping unused ones)
• For each numeric predictor (covariate), use its average.¹

The reference grid is then a regular grid of all combinations of these reference levels.

As a simple example, consider again the pigs dataset (see help("fiber") for details). Examination of residual plots from preliminary models suggests that it is a good idea to work in terms of log concentration.

If we treat the predictor percent as a factor, we might fit the following model:

pigs.lm1 <- lm(log(conc) ~ source + factor(percent), data = pigs)

The reference grid for this model can be found via the ref_grid function:

ref_grid(pigs.lm1)

## 'emmGrid' object with variables:
##     source = fish, soy, skim
##     percent =  9, 12, 15, 18
## Transformation: "log"

(Note: Many of the calculations that follow are meant to illustrate what is inside this reference-grid object; You don’t need to do such calculations yourself in routine analysis; just use the emmeans() (or possibly ref_grid()) function as we do later.)

In this model, both predictors are factors, and the reference grid consists of the \(3\times4 = 12\) combinations of these factor levels. It can be seen explicitly by looking at the grid slot of this object:

ref_grid(pigs.lm1) @ grid

##    source percent .wgt.
## 1    fish       9     2
## 2     soy       9     3
## 3    skim       9     3
## 4    fish      12     3
## 5     soy      12     3
## 6    skim      12     3
## 7    fish      15     2
## 8     soy      15     3
## 9    skim      15     2
## 10   fish      18     3
## 11    soy      18     1
## 12   skim      18     1

Note that other information is retained in the reference grid, e.g., the transformation used on the response, and the cell counts as the .wgt. column.

Now, suppose instead that we treat percent as a numeric predictor. This leads to a different model – and a different reference grid.

pigs.lm2 <- lm(log(conc) ~ source + percent, data = pigs)
ref_grid(pigs.lm2)

## 'emmGrid' object with variables:
##     source = fish, soy, skim
##     percent = 12.931
## Transformation: "log"

This reference grid has the levels of source, but only one percent value, its average. Thus, the grid has only three elements:

ref_grid(pigs.lm2) @ grid

##   source  percent .wgt.
## 1   fish 12.93103    10
## 2    soy 12.93103    10
## 3   skim 12.93103     9

Back to Contents

Estimated marginal means

Once the reference grid is established, we can consider using the model to estimate the mean at each point in the reference grid. (Curiously, the convention is to call this “prediction” rather than “estimation”). For pigs.lm1, we have

pigs.pred1 <- matrix(predict(ref_grid(pigs.lm1)), nrow = 3)
pigs.pred1

##          [,1]     [,2]     [,3]     [,4]
## [1,] 3.220292 3.399846 3.437691 3.520141
## [2,] 3.493060 3.672614 3.710459 3.792909
## [3,] 3.622569 3.802124 3.839968 3.922419

Estimated marginal means (EMMs) are defined as equally weighted means of these predictions at specified margins:

apply(pigs.pred1, 1, mean) ### EMMs for source

## [1] 3.394492 3.667260 3.796770

apply(pigs.pred1, 2, mean) ### EMMs for percent

## [1] 3.445307 3.624861 3.662706 3.745156

For the other model, pigs.lm2, we have only one point in the reference grid for each source level; so the EMMs for source are just the predictions themselves:

predict(ref_grid(pigs.lm2))

## [1] 3.379865 3.652693 3.783120

These are slightly different from the previous EMMs for source, emphasizing the fact that EMMs are model-dependent. In models with covariates, EMMs are often called adjusted means.

The emmeans function computes EMMs, accompanied by standard errors and confidence intervals. For example,

emmeans(pigs.lm1, "percent")

##  percent emmean     SE df lower.CL upper.CL
##        9   3.45 0.0409 23     3.36     3.53
##       12   3.62 0.0384 23     3.55     3.70
##       15   3.66 0.0437 23     3.57     3.75
##       18   3.75 0.0530 23     3.64     3.85
## 
## Results are averaged over the levels of: source 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

In these examples, all the results are presented on the log(conc) scale (and the annotations in the output warn of this). It is possible to convert them back to the conc scale by back-transforming. This topic is discussed in the vignette on transformations.

An additional note: There is an exception to the definition of EMMs given here. If the model has a nested structure in the fixed effects, then averaging is performed separately in each nesting group. See the section on nesting in the “messy-data” vignette for an example.

Back to Contents

Altering the reference grid

It is possible to alter the reference grid. We might, for example, want to define a reference grid for pigs.lm2 that is comparable to the one for pigs.lm1.

ref_grid(pigs.lm2, cov.keep = "percent")

## 'emmGrid' object with variables:
##     source = fish, soy, skim
##     percent =  9, 12, 15, 18
## Transformation: "log"

Using cov.keep = "percent" specifies that, instead of using the mean, the reference grid should use all the unique values of each covariate“percent”`.

Another option is to specify a cov.reduce function that is used in place of the mean; e.g.,

ref_grid(pigs.lm2, cov.reduce = range)

## 'emmGrid' object with variables:
##     source = fish, soy, skim
##     percent =  9, 18
## Transformation: "log"

Another option is to use the at argument. Consider this model for the built-in mtcars dataset:

mtcars.lm <- lm(mpg ~ disp * cyl, data = mtcars)
ref_grid(mtcars.lm)

## 'emmGrid' object with variables:
##     disp = 230.72
##     cyl = 6.1875

Since both predictors are numeric, the default reference grid has only one point. For purposes of describing the fitted model, you might want to obtain predictions at a grid of points, like this:

mtcars.rg <- ref_grid(mtcars.lm, cov.keep = 3,
                      at = list(disp = c(100, 200, 300)))
mtcars.rg

## 'emmGrid' object with variables:
##     disp = 100, 200, 300
##     cyl = 4, 6, 8

This illustrates two things: a new use of cov.keep and the at argument. cov.keep = "3" specifies that any covariates having 3 or fewer unique values is treated like a factor (the system default is cov.keep = "2"). The at specification gives three values of disp, overriding the default behavior to use the mean of disp. Another use of at is to focus on only some of the levels of a factor. Note that at does not need to specify every predictor; those not mentioned in at are handled by cov.reduce, cov.keep, or the default methods. Also, covariate values in at need not be values that actually occur in the data, whereas cov.keep will use only values that are achieved.

Back to Contents

Derived covariates

You need to be careful when one covariate depends on the value of another. To illustrate in the mtcars example, suppose we want to use cyl as a factor and include a quadratic term for disp:

mtcars.1 <- lm(mpg ~ factor(cyl) + disp + I(disp^2), data = mtcars)
emmeans(mtcars.1, "cyl")

##  cyl emmean   SE df lower.CL upper.CL
##    4   19.3 2.66 27     13.9     24.8
##    6   17.2 1.36 27     14.4     20.0
##    8   18.8 1.47 27     15.7     21.8
## 
## Confidence level used: 0.95

Some users may not like function calls in the model formula, so they instead do something like this:

mtcars <- transform(mtcars, 
                    Cyl = factor(cyl),
                    dispsq = disp^2)
mtcars.2 <- lm(mpg ~ Cyl + disp + dispsq, data = mtcars)
emmeans(mtcars.2, "Cyl")

##  Cyl emmean   SE df lower.CL upper.CL
##  4     20.8 2.05 27     16.6     25.0
##  6     18.7 1.19 27     16.3     21.1
##  8     20.2 1.77 27     16.6     23.9
## 
## Confidence level used: 0.95

Wow! Those are really different results – even though the models are equivalent. Why is this? To understand, look at the reference grids:

ref_grid(mtcars.1)

## 'emmGrid' object with variables:
##     cyl = 4, 6, 8
##     disp = 230.72

ref_grid(mtcars.2)

## 'emmGrid' object with variables:
##     Cyl = 4, 6, 8
##     disp = 230.72
##     dispsq = 68113

For both models, the reference grid uses the disp mean of 230.72. But for mtcars.2, we also set dispsq to its mean of 68113. This is not right, because dispsq should be the square of disp (about 53232, not 68113) in order to be consistent. If we use that value of dispsq, we get the same results (modulus rounding error) as for mtcars.1:

emmeans(mtcars.2, "Cyl", at = list(dispsq = 230.72^2))

##  Cyl emmean   SE df lower.CL upper.CL
##  4     19.3 2.66 27     13.9     24.8
##  6     17.2 1.36 27     14.4     20.0
##  8     18.8 1.47 27     15.7     21.8
## 
## Confidence level used: 0.95

In summary, for polynomial models and others where some covariates depend on others in nonlinear ways, include that dependence in the model formula (as in mtcars.1) using I() or poly() expressions, or alter the reference grid so that the dependency among covariates is correct.

Non-predictor variables

Reference grids are derived using the variables in the right-hand side of the model formula. But sometimes, these variables are not actually predictors. For example:

deg <- 2
mod <- lm(y ~ treat * poly(x, degree = deg), data = mydata)

If we call ref_grid() or emmeans() with this model, it will try to construct a grid of values of treat, x, and deg – causing an error because deg is not a predictor in this model. To get things to work correctly, you need to name deg in a params argument, e.g.,

emmeans(mod, ~ treat | x, at = list(x = 1:3), params = "deg")

Back to Contents

Graphical displays

The results of ref_grid() or emmeans() (these are objects of class emmGrid) may be plotted in two different ways. One is an interaction-style plot, using emmip(). In the following, let’s use it to compare the predictions from pigs.lm1 and pigs.lm2:

emmip(pigs.lm1, source ~ percent)

emmip(ref_grid(pigs.lm2, cov.reduce = FALSE), source ~ percent)

Notice that emmip() may also be used on a fitted model. The formula specification needs the x variable on the right-hand side and the “trace” factor (what is used to define the different curves) on the left. This is a good time to yet again emphasize that EMMs are based on a model. Neither of these plots is an interaction plot of the data; they are interaction plots of model predictions; and since both models do not include an interaction, no interaction at all is evident in the plots.

plot(mtcars.rg, by = "disp")

mtcars.rg_d.c <- ref_grid(mtcars.lm, at = list(cyl = c(4,6,8)),
                          cov.reduce = disp ~ cyl)
mtcars.rg_d.c @ grid

##        disp cyl .wgt.
## 1  93.78673   4     1
## 2 218.98458   6     1
## 3 344.18243   8     1

plot(mtcars.rg_d.c)

require("ggplot2")

## Loading required package: ggplot2

emmip(pigs.lm1, ~ percent | source, CIs = TRUE) +
    geom_point(aes(x = percent, y = log(conc)), data = pigs, pch = 2, color = "blue")

Formatting results

If you want to include emmeans() results in a report, you might want to have it in a nicer format than just the printed output. We provide a little bit of help for this, especially if you are using RMarkdown or SWeave to prepare the report. There is an xtable method for exporting these results, which we do not illustrate here but it works similarly to xtable() in other contexts. Also, the export option the print() method allows the user to save exactly what is seen in the printed output as text, to be saved or formatted as the user likes (see the documentation for print.emmGrid for details). Here is an example using one of the objects above:

ci <- confint(mtcars.rg_d.c, level = 0.90, adjust = "scheffe")
xport <- print(ci, export = TRUE)
cat("<font color = 'blue'>\n")
knitr::kable(xport$summary, align = "r")
for (a in xport$annotations) cat(paste(a, "<br>"))
cat("</font>\n")

Confidence level used: 0.9
Conf-level adjustment: scheffe method with rank 3

Back to Contents

Using weights

It is possible to override the equal-weighting method for computing EMMs. Using weights = "cells" in the call will weight the predictions according to their cell frequencies (recall this information is retained in the reference grid). This produces results comparable to ordinary marginal means:

emmeans(pigs.lm1, "percent", weights = "cells")

##  percent emmean     SE df lower.CL upper.CL
##        9   3.47 0.0407 23     3.39     3.56
##       12   3.62 0.0384 23     3.55     3.70
##       15   3.67 0.0435 23     3.58     3.76
##       18   3.66 0.0515 23     3.55     3.76
## 
## Results are averaged over the levels of: source 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

Note that, as in the ordinary means in the motivating example, the highest estimate is for percent = 15 rather than percent = 18. It is interesting to compare this with the results for a model that includes only percent as a predictor.

pigs.lm3 <- lm(log(conc) ~ factor(percent), data = pigs)
emmeans(pigs.lm3, "percent")

##  percent emmean     SE df lower.CL upper.CL
##        9   3.47 0.0731 25     3.32     3.62
##       12   3.62 0.0689 25     3.48     3.77
##       15   3.67 0.0782 25     3.51     3.83
##       18   3.66 0.0925 25     3.46     3.85
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

The EMMs in these two tables are identical, but their standard errors are considerably different. That is because the model pigs.lm1 accounts for variations due to source. The lesson here is that it is possible to obtain statistics comparable to ordinary marginal means, while still accounting for variations due to the factors that are being averaged over.

Back to Contents

Multivariate responses

The emmeans package supports various multivariate models. When there is a multivariate response, the dimensions of that response are treated as if they were levels of a factor. For example, the MOats dataset provided in the package has predictors Block and Variety, and a four-dimensional response yield giving yields observed with varying amounts of nitrogen added to the soil. Here is a model and reference grid:

MOats.lm <- lm (yield ~ Block + Variety, data = MOats)
ref_grid (MOats.lm, mult.name = "nitro")

## 'emmGrid' object with variables:
##     Block = VI, V, III, IV, II, I
##     Variety = Golden Rain, Marvellous, Victory
##     nitro = multivariate response levels: 0, 0.2, 0.4, 0.6

So, nitro is regarded as a factor having 4 levels corresponding to the 4 dimensions of yield. We can subsequently obtain EMMs for any of the factors Block, Variety, nitro, or combinations thereof. The argument mult.name = "nitro" is optional; if it had been excluded, the multivariate levels would have been named rep.meas.

Back to Contents

A set of comparisons or well-chosen contrasts is more useful and interpretable than an omnibus F test

F tests are useful for model selection, but don’t tell you anything specific about the nature of an effect. If F has a small P value, it suggests that there is some effect, somewhere. It doesn’t even necessarily imply that any two means differ statistically.

Use adjusted P values

When you run a bunch of tests, there is a risk of making too many type-I errors, and adjusted P values (e.g., the Tukey adjustment for pairwise comparisons) keep you from making too many mistakes. That said, it is possible to go overboard; and it’s usually reasonable to regard each “by” group as a separate family of tests for purposes of adjustment.

It is not necessary to have a significant F test as a prerequisite to doing comparisons or contrasts

… as long as an appropriate adjustment is used. There do exist rules such as the “protected LSD” by which one is given license to do unadjusted comparisons provided the \(F\) statistic is “significant.” However, this is a very weak form of protection for which the justification is, basically, “if \(F\) is significant, you can say absolutely anything you want.”

Get the model right first

Everything the emmeans package does is an interpretation of the model that you fitted to the data. If the model is bad, you will get bad results from emmeans() and other functions. Every single limitation of your model, be it presuming constant error variance, omitting interaction terms, etc., becomes a limitation of the results emmeans() produces. So do a responsible job of fitting the model. And if you don’t know what’s meant by that…

Consider seeking the advice of a statistical consultant

Statistics is hard. It is a lot more than just running programs and copying output. It is your research; is it important that it be done right? Many academic statistics and biostatistics departments can refer you to someone who can help.

Back to Contents

References

Wasserman RL, Lazar NA (2016) “The ASA’s Statement on p-Values: Context, Process, and Purpose,” The American Statistician, 70, 129–133, https://doi.org/10.1080/00031305.2016.1154108

Wasserman RL, Schirm AL, Lazar, NA (2019) “Moving to a World Beyond ‘p < 0.05’,” The American Statistician, 73, 1–19, https://doi.org/10.1080/00031305.2019.1583913