FAQs for emmeans

emmeans package, Version 1.4.1

This vignette contains answers to questions received from users or posted on discussion boards like Cross Validated and Stack Overflow

Contents

  1. What are EMMs/lsmeans?
  2. What is the fastest way to obtain EMMs and pairwise comparisons?
  3. I have three (or two or four) factors that interact
  4. I have covariate(s) that interact(s) with factor(s)
  5. I have covariate(s) and am fitting a polynomial model
  6. Some “significant” comparisons have overlapping confidence intervals
  7. All my pairwise comparisons have the same P value
  8. emmeans() doesn’t work as expected
  9. All or some of the results are NA
  10. If I analyze subsets of the data separately, I get different results
  11. My lsmeans/EMMs are way off from what I expected
  12. Why do I get Inf for the degrees of freedom?
  13. I get exactly the same comparisons for each “by” group
  14. My ANOVA F is significant, but no pairwise comparisons are
  15. I asked for a Tukey adjustments, but that’s not what I got
  16. emmeans() gives me pooled t tests, but I expected Welch’s t

Index of all vignette topics

What are EMMs/lsmeans?

Estimated marginal means (EMMs), a.k.a. least-squares means, are predictions on a reference grid of predictor settings, or marginal averages thereof. See details in the “basics” vignette.

What is the fastest way to obtain EMMs and pairwise comparisons?

There are two answers to this (i.e., be careful what you wish for):

  1. Don’t think; just fit the first model that comes to mind and run emmeans(model, pairwise ~ treatment). This is the fastest way; however, the results have a good chance of being invalid.
  2. Do think: Make sure you fit a model that really explains the responses. Do diagnostic residual plots, include appropriate interactions, account for heteroscadesticity if necessary, etc. This is the fastest way to obtain appropriate estimates and comparisons.

The point here is that emmeans() summarizes the model, not the data directly. If you use a bad model, you will get bad results. And if you use a good model, you will get appropriate results. It’s up to you: it’s your research—is it important?

Back to Contents

I have three (or two or four) factors that interact

Perhaps your question has to do with interacting factors, and you want to do some kind of post hoc analysis comparing levels of one (or more) of the factors on the response. Some specific versions of this question…

My first answer is: plots almost always help. If you have factors A, B, and C, try something like emmip(model, A ~ B | C), which creates an interaction-style plot of the predictions against B, for each A, with separate panels for each C. This will help visualize what effects stand out in a practical way. This can guide you in what post-hoc tests would make sense. See the “interactions” vignette for more discussion and examples.

Back to Contents

I have covariate(s) and am fitting a polynomial model

You need to be careful to define the reference grid consistently. For example, if you use covariates x and xsq (equal to x^2) to fit a quadratic curve, the default reference grid uses the mean of each covariate – and mean(xsq) is usually not the same as mean(x)^2. So you need to use at to ensure that the covariates are set consistently with respect to the model. See this subsection of the “basics” vignette for an example.

3. Some “significant” comparisons have overlapping confidence intervals

That can happen because it is just plain wrong to use [non-]overlapping CIs for individual means to do comparisons. Look at the printed results from something like emmeans(mymodel, pairwise ~ treatment). In particular, note that the SE values are not the same*, and may even have different degrees of freedom. Means are one thing statistically, and differences of means are quite another thing. Don’t ever mix them up, and don’t ever use a CI display for comparing means.

I’ll add that making hard-line decisions about “significant” and “non-significant” is in itself a poor practice. See the discussion in the “basics” vignette

All my pairwise comparisons have the same P value

This will happen if you fitted a model where the treatments you want to compare were put in as a numeric predictor; for example dose, with values of 1, 2, and 3. If dose is modeled as numeric, you will be fitting a linear trend in those dose values, rather than a model that allows those doses to differ in arbitrary ways. Go back and fit a different model using factor(dose) instead; it will make all the difference. This is closely related to the next topic.

emmeans() doesn’t work as expected

Equivalently, users ask how to get post hoc comparisons when we have covariates rather than factors. Yes, it does work, but you have to tell it the appropriate reference grid.

But before saying more, I have a question for you: Are you sure your model is meaningful?

Assuming that issue is settled, you can do something like emmeans(model, "sex", at = list(sex = 1:2)), to get get separate EMMs for each sex rather than one EMM for the average numeric value of sex.

An alternative to the at list is to use cov.reduce = FALSE, which specifies that the unique values of covariates are to be used rather than reducing them to their means. However, the specification applies to all covariates, so if you have another one, say age, that has 43 different values in your data, you will have a mess on your hands.

See “altering the reference grid” in the “basics” vignette for more discussion.

Back to Contents

All or some of the results are NA

The emmeans package uses tools in the estimability package to determine whether its results are uniquely estimable. For example, in a two-way model with interactions included, if there are no observations in a particular cell (factor combination), then we cannot estimate the mean of that cell.

When some of the EMMs are estimable and others are not, that is information about missing information in the data. If it’s possible to remove some terms from the model (particularly interactions), that may make more things estimable if you re-fit with those terms excluded; but don’t delete terms that are really needed for the model to fit well.

When all of the estimates are non-estimable, it could be symptomatic of something else. Some possibilities include:

pg <- transform(pigs, x = rep(1:3, c(10, 10, 9)))
pg.lm <- lm(log(conc) ~ x + source + factor(percent), data = pg)
emmeans(pg.lm, consec ~ percent)
## $emmeans
##  percent emmean SE df asymp.LCL asymp.UCL
##        9 nonEst NA NA        NA        NA
##       12 nonEst NA NA        NA        NA
##       15 nonEst NA NA        NA        NA
##       18 nonEst NA NA        NA        NA
## 
## Results are averaged over the levels of: source 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate     SE df t.ratio p.value
##  12 - 9     0.1796 0.0561 23 3.202   0.0108 
##  15 - 12    0.0378 0.0582 23 0.650   0.8613 
##  18 - 15    0.0825 0.0691 23 1.194   0.5200 
## 
## Results are averaged over the levels of: source 
## Results are given on the log (not the response) scale. 
## P value adjustment: mvt method for 3 tests

The “messy-data” vignette has more examples and discussion.

Back to Contents

If I analyze subsets of the data separately, I get different results

Estimated marginal means summarize the model that you fitted to the data – not the data themselves. Many of the most common models rely on several simplifying assumptions – that certain effects are linear, that the error variance is constant, etc. – and those assumptions are passed forward into the emmeans() results. Doing separate analyses on subsets usually comprises departing from that overall model, so of course the results are different.

My lsmeans/EMMs are way off from what I expected

First step: Carefully read the annotations below the output. Do they say something like “results are on the log scale, not the response scale”? If so, that explains it. A Poisson or logistic model involves a link function, and by default, emmeans() produces its results on that same scale. You can add type = "response" to the emmeans() call and it will put the results of the scale you expect. But that is not always the best approach. The “transformations” vignette has examples and discussion.

Why do I get Inf for the degrees of freedom?

This is simply the way that emmeans labels asymptotic results (that is, estimates that are tested against the standard normal distribution – z tests – rather than the t distribution). Note that obtaining quantiles or probabilities from the t distribution with infinite degrees of freedom is the same as obtaining the corresponding values from the standard normal. For example:

qt(c(.9, .95, .975), df = Inf)
## [1] 1.281552 1.644854 1.959964
qnorm(c(.9, .95, .975))
## [1] 1.281552 1.644854 1.959964

so when you see infinite d.f., that just means its a z test or a z confidence interval.

Back to Contents

I get exactly the same comparisons for each “by” group

As mentioned elsewhere, EMMs summarize a model, not the data. If your model does not include any interactions between the by variables and the factors for which you want EMMs, then by definition, the effects for the latter will be exactly the same regardless of the by variable settings. So of course the comparisons will all be the same. If you think they should be different, then you are saying that your model should include interactions between the factors of interest and the by factors.

My ANOVA F is significant, but no pairwise comparisons are

First of all, you should not be making binary decisions of “significant” or “nonsignificant.” This is a simplistic view of P values that assigns an unmerited magical quality to the value 0.05. It is suggested that you just report the P values actually obtained, and let your readers decide how significant your findings are in the context of the scientific findings.

But to answer the question: This is a common misunderstanding of ANOVA. If F has a particular P value, this implies only that some contrast among the means (or effects) has the same P value, after applying the Scheffe adjustment. That contrast may be very much unlike a pairwise comparison, especially when there are several means being compared. Having an F statistic with a P value of, say, 0.06, does not imply that any pairwise comparison will have a P value of 0.06 or smaller. Again referring to the paragraph above, just report the P value for each pairwise comparison, and don’t try to relate them to the F statistic.

Another consideration is that by default, P values for pairwise comparisons are adjusted using the Tukey method, and the adjusted P values can be quite a bit larger than the unadjusted ones. (But I definitely do not advocate using no adjustment to “repair” this problem.)

I asked for Tukey adjustments, but that’s not what I got

There are two reasons this could happen:

  1. There are only two means (or only two in each by group). Thus there is only one comparison. When there is only one thing to test, there is no multiplicity issue, and hence no multiplicity adjustment to the P values.
  2. A Tukey adjustment is inappropriate. The Tukey adjustment is appropriate for pairwise comparisons of means. When you have some other set of contrasts, the Tukey method is deemed unsuitable and the Sidak method is used instead. A suggestion is to use "mvt" adjustment (which is exact); we don’t default to this because it can require a lot of computing time for a large set of contrasts or comparisons.

Back to Contents

emmeans() gives me pooled t tests, but I expected Welch’s t

It is important to note that emmeans() and its relatives produce results based on the model object that you provide – not the data. So if your sample SDs are wildly different, a model fitted using lm() or aov() is not a good model, because those R functions use a statistical model that presumes that the errors have constant variance. That is, the problem isn’t in emmeans(), it’s in handing it an inadequate model object.

Here is a simple illustrative example. Consider a simple one-way experiment and the following model:

mod1 <- aov(response ~ treat, data = mydata)
emmeans(mod1, pairwise ~ treat)

This code will estimate means and comparisons among treatments. All standard errors, confidence intervals, and t statistics are based on the pooled residual SD with N - k degrees of freedom (assuming N observations and k treatments). These results are useful only if the underlying assumptions of mod1 are correct – including the assumption that the error SD is the same for all treatments.

Alternatively, you could fit the following model using generalized least-squares:

mod2 = nlme::gls(response ~ treat, data = mydata,
                 weights = varIdent(form = ~1 | treat))
emmeans(mod2, pairwise ~ treat)

This model specifies that the error variance depends on the levels of treat. This would be a much better model to use when you have wildly different sample SDs. The results of the emmeans() call will reflect this improvement in the modeling. The standard errors of the EMMs will depend on the individual sample variances, and the t tests of the comparisons will be in essence the Welch t statistics with Satterthwaite degrees of freedom.

To obtain appropriate post hoc estimates, contrasts, and comparisons, one must first find a model that successfully explains the peculiarities in the data. This point cannot be emphasized enough. If you give emmeans() a good model, you will obtain correct results; if you give it a bad model, you will obtain incorrect results. Get the model right first.

Back to Contents

Index of all vignette topics