0

Here is my data frame (my real DF has way more data points):

    rearing.temp<-c("15", "15", "15", "15", "19", "19", "19", "19")
    source<-c("field", "field", "woods", "woods", "field", "field", "woods", "woods")
    runway.temp<-c("40","20","40","20","40","20","40","20")
    velocity<-c("2.3", "2.1", "1.9", "1.9", "2.3", "2.2", "2.3", "2.0")
    snail<-data.frame(rearing.temp, source, runway.temp, velocity)

Here is my model:

mod <- lmer(velocity ~ runway.temp*source*rearing.temp + (1|family) + (1|collection site) + (1|individual.plus.family.id.combined), data=snail)

When I do an emmeans contrast:

emmeans(mod, pairwise~runway.temp*source*rearing.temp)

I get 28 different comparisons, but I am only interested in looking at the difference between the velocity of field snails reared at 15° tested at the 40° runway temperature compared to woods snails reared at 15° tested at the 40° runway temperature. I just want to do one comparison between the snails that are both reared at the same temp, are tested at the same temp, but are sourced from different habitats. How can I do this?

Thank you,

Ricardo

Ricardo Yomez
  • 11
  • 3
  • Why are you using a model when you only have one observation per group? Doing emmeans is just going to give you the difference between pairs of values. Can you not calculate the differences directly from the data frame? – Kevin A Dec 16 '20 at 19:58
  • I should have clarified: My real dataset has ~800 observations, I just wanted to make something reproducible to get my idea across – Ricardo Yomez Dec 17 '20 at 00:44
  • Anyone? I could really use some help – Ricardo Yomez Dec 20 '20 at 01:41
  • You are fitting a simple linear model with a three-way interaction term to your data. This has nothing to do with nested levels. If you have a nested experimental design, you may want to fit a nested/mixed effect model (e.g. using `lme4` or `lmer`). Regarding pairwise comparisons: You must perform *all* pairwise comparisons and not just the one that interests you. – Maurits Evers Jan 17 '21 at 23:50
  • I recommend consulting with a resident statistician and have him explain to your boss why manually selecting comparisons is a dangerous game to play. This is the kind of selective p-value fishing that has come under a lot of criticism. Your boss needs to understand that, yes, you absolutely need to "*compare every single level against the other*" and, yes, that will be *"affecting the p-values"*. This is exactly the point of multiple comparisons in post-hoc tests. It is wrong (i.e. statistically unsound) to *[d]o planned contrasts of only the levels you are interested in comparing.* – Maurits Evers Jan 18 '21 at 22:29
  • [continued] You (and your boss) can find many posts discussing this on Stack Exchange (specifically on Cross Validated, the statistics platform); there are also many papers in academic journals discussing this form of statistical malpractice in various research areas. The bottom line is: Discussing this on an internet forum is probably not the right place. Your best course of action in my opinion is to engage a local statistician or statistical consultant and help him to help you. – Maurits Evers Jan 18 '21 at 22:33
  • This guy seems to think it's ok https://aosmith.rbind.io/2019/04/15/custom-contrasts-emmeans/ – Ricardo Yomez Jan 31 '21 at 00:54
  • Maurits, I don't think you understand my motives for doing a selective comparison: Comparing snails run at 40 degrees to snails run at 20 degrees is meaningless. I don't need that contrast because the contrast is only meaningful for snails within the same runway temperature and rearing temperature treatment. The whole point of my experiment is to look at differences between source habitats. – Ricardo Yomez Jan 31 '21 at 01:16

1 Answers1

1

The contrast function provides for specifying a named list of desired contrast coefficients. In your case, apparently you have 8 means (2x2x2 array), and they are arranged in order with the first index alternating fastest and the last factor the slowest - (1,1,1), (2,1,1), (1,2,1), ..., (2,2,2). So do something like this:

emm <- emmeans(mod, ~runway.temp*source*rearing.temp)
custom <- list(`111vs211` = c(1,0,0,0,-1,0,0,0),
               `211vs121` = c(0,1,-1,0,0,0,0,0))
contrast(emm, custom)

There is abundant documentation and examples. See for example ? contrast and vignette("contrasts", "emmeans")

Russ Lenth
  • 5,922
  • 2
  • 13
  • 21
  • Thank you for your help! I get the error, "Error: unexpected symbol in "custom <- list(111vs211" when I try to run your code. If I replace 111vs211 (which I think is the problem) with A2 and B2, I get the error, "Error in UseMethod("contrast") : no applicable method for 'contrast' applied to an object of class "function." I think if I can run this code, you've solved the problem! – Ricardo Yomez Jan 31 '21 at 01:01
  • That means you didn't create emm first. Do all the steps. For the first issue, you need to add the back-quotes shown. I did fail to include those at first – Russ Lenth Jan 31 '21 at 02:34
  • Thanks Russ for your response. I'm very close; I now get: Error in contrast.emmGrid(emm, custom) : Nonconforming number of contrast coefficients. I think I just need to work on understanding the concatenated part. What do you mean when you say I have "8 means"? I will need to read up on what a 2x2x2 array is. I'm unclear on how my variables get translated into this array, and how parts of that array can be pulled out by using the concatenate function & using either 1, 0, or -1. I tried doing 'vignette("contrasts", "emmeans")' but in V4.0.2 it doesn't work. Read ? contrast though still unclear – Ricardo Yomez Jan 31 '21 at 14:43
  • My real data has 66 contrasts, so wouldn't I need 66 zeros inside the concatenate with whichever ones I want as 1s? – Ricardo Yomez Jan 31 '21 at 14:55
  • How many estimates are displayed when you list `emm`? You need exactly that many coefficients in each set. – Russ Lenth Jan 31 '21 at 18:49
  • 3
    I think there may be a serious disconnect here in terms of concepts. I suggest finding consulting help, perhaps from a nearby university statistics department. Statistics is a whole lot more than getting programs to run. – Russ Lenth Jan 31 '21 at 18:51