3

I hope somebody is available to help a desperate rookie..

I fitted a glmer with a Poisson distribution and log link, including main effects and several interactions, an offset variable and a random effect. Something like this:

model1 <- glmer(count ~ A + B + C + A*B + offset(log(X)) + ( 1 | D), data = workmap, family = poisson(link = "log"))

A, B and C are all factors with two levels.

My Prof generated "contrasts of marginal linear predictions" in Stata to, for example, look at the contrasts provided by A@B, or just simply A. I now want to do the same but in R by making use of the emmeans package.

First: should I use emmeans() or contrast() command? What is the difference? I did the following: emm <- emmeans(model1, pairwise ~ A | B)

Was that the right thing to do or should I have used the contrast() command?

Ronja
  • 41
  • 4

2 Answers2

3

The emmeans package provides some flexibility in looking at different parts of the analysis, as well as some convenience functions. With this example, you could do:

EMM <- emmeans(model1, ~ A|B)
EMM          # show the estimated marginal means

PRS <- contrast(EMM, "pairwise")  # or PRS <- pairs(EMM)
PRS          # show the pairwise comparisons

The call you have, emmeans(model1, pairwise ~ A|B), does both of the above, and packages the results in a list. So in that sense, there is no difference.

That said, I strongly recommend doing the estimation of EMMs, and contrasts thereof, separately, because there may be more things you want to do with either part. For example,

# show on the original scale rather than the log scale:
confint(EMM, type = "response")

# show estimated Poisson rates:
confint(EMM, offset = 0)

# show pairwise comparisons as ratios on original scale:
summary(PRS, type = "response", infer = c(TRUE, TRUE))

Also, the fact that the emmeans(..., pairwise ~ ...) construct creates a list of emmGrid objects rather than a single emmGrid object causes confusion for some users.

By the way, since you have a mixed model, there is an additional issue that back-transformed estimates (with type = "response") are biased. To correct for bias, you need an estimate of the SD of the random effect, which you can get from VarCorr(model1). Then you can do confint(EMM, type = "response", bias.adj = TRUE, sigma = <your SD estimate>). It's important to specify sigma because the default value is incorrect for these models.

Russ Lenth
  • 15,161
  • 20
  • 53
  • wow, you already helped me a lot with your answer. Your step by step description is fantastic. I hope, I may ask you a follow-up question: I am mainly interested in the pairwise comparison (PRS). Can I run the confint command with PRS instead of EMM? – Ronja Oct 01 '21 at 18:10
  • Sorry to flood you with another question... which output should I use to report contrasts (contrast(EMM/PRS, "pairwise"), confint(EMM, type = "response") once corrected or confint(EMM/PRS, offset = 0)). I would like to be able to say something like "after implementation of treatment (A), dogs (B) reduced their barking by 35% compared to wolves (B) which only reduced barking by 2%. – Ronja Oct 01 '21 at 18:21
  • confint() or test() with either object return CIs or tests, response. summary(..., infer = c(TRUE, TRUE)) shows both. – Russ Lenth Oct 01 '21 at 18:58
  • Other Q depends on context. With no offset, you're including X as a mediating variable, which likely has a different mean with each treatment. With offset = 0 (= log(1)), you're comparing barking assuming the same X = 1 for each treatment. The latter is probably more meaningful unless you think the treatment causes changes in X as well as in barking. – Russ Lenth Oct 01 '21 at 19:04
  • Thank you very much, I guess "ratio" is the Odds ration? Is there any possibility to replace it with the Incidence rate ratio (IRR)? – Ronja Oct 01 '21 at 19:33
  • It is the IRR already. We can't get odds ratios except with a binomial logit model. Again, you should consider whether or not you want offset = 0 or not in those ratios. – Russ Lenth Oct 01 '21 at 20:09
  • I figured it out, thanks to you! Only problem left: A has to levels (0 and 1) and in regard to the contrasts, I get the ratio for level 0 relative to level 1, but would like the inverse. Any advice? – Ronja Oct 05 '21 at 16:43
  • Try "revpairwise" – Russ Lenth Oct 05 '21 at 19:57
  • works like a charm! Thank you, from the bottom of my heart. – Ronja Oct 05 '21 at 21:29
2

Classically in statistics

a contrast is a linear combination of variables (parameters or statistics) whose coefficients add up to zero, allowing comparison of different treatments.

In practice, the term has often been expanded to include any linear combination whether or not the coefficients of the contrast add up to zero, as noted in the emmeans vignette on contrasts.

If you have the coefficient estimates and their covariance matrix, you can use the formulas for the variance of a linear combination to calculate any desired contrast (in that extended sense) and its standard error yourself. If you're just learning, it might make sense to try that yourself on some simple cases.

The emmeans package provides some very useful, widely applicable tools to make such calculations easier. It does that by first constructing what's called a "reference grid" for the underlying model, as explained in a vignette. The model is based on the data directly while the reference grid is based on the model of the data.

So the first step in using any of the emmeans tools is to set up the reference grid. That's done conveniently as part of a call to the emmeans() function. The way you called the function, however:

emm <- emmeans(model1, pairwise ~ A | B)

also generated the indicated set of pairwise contrasts based on the reference grid it generated. If you look at your emm object you will find that it is a list of 2 emmGrid objects, one for the specified aspects of the original model and a second representing the specified contrasts.

In emmeans the contrast() function only works on an emmGrid object. That allows you to evaluate additional contrasts beyond what you first considered in your call to emmeans() without having to rebuild the grid for the original model. It also allows you to work with derived grids. For example, yesterday I did further contrasts on a grid representing an initial set of contrasts.

That's not the way that a contrast() function works in all software, however. In the rms package the contrast() function works directly on models. And those type of contrast() functions are different from the contrasts() function in the standard R stats package that reports or specifies how a particular categorical variable will be coded in a model in the first place.

So read very carefully to see how the word "contrast" and any associated functions are being used in a given context.

EdM
  • 57,766
  • 7
  • 66
  • 187