2
$\begingroup$

I am trying to analyse subgroups of a mixed-effects model using custom contrasts. I have 2 factors F1 and F2, each with 3 levels. My interest is in comparing the groups like this:

  1. F1_group1 vs. F1_group2and3
  2. F1_group2 vs. F1_group3
  3. F2_group1 vs. F2_group2and3
  4. F2_group2 vs. F2_group3

I am using lme4::lmer(), and them emmeans::emmeans(). The model looks good, the contrasts look good (uncorrelated etc.). The things that "should" be significant are, and those that "should not" are not.

Problem: The estimates are obviously scaled. The exact values are way too large. If I do paired comparisons the estimates are fine. I do not understand how exactly and can't find information on this anywhere.

I am very interested in the exact estimates, not just the p-values, so I would like to figure out where the issue occurs. Any ideas or pointers?

Contrasts look like this:

F1 F2 c1 c2 c3 c4
a a -2 0 1 -1
a b -2 0 -2 0
a c -2 0 1 1
b a 1 -1 1 -1
b b 1 -1 -2 0
b c 1 -1 1 1
c a 1 1 1 -1
c b 1 1 -2 0
c c 1 1 1 1
$\endgroup$
3
  • $\begingroup$ How do you know the estimates are too large ? Can you include the output that shows this please ? $\endgroup$ Commented Jun 5, 2021 at 8:52
  • $\begingroup$ Thank you so much for the question, it got me on the right track. Some scaling was indeed required. $\endgroup$
    – c06n
    Commented Jun 7, 2021 at 11:35
  • $\begingroup$ You're welcome. Glad I could help ! $\endgroup$ Commented Jun 7, 2021 at 11:36

2 Answers 2

2
$\begingroup$

My guess is that your expectation is that the regression coefficients will be estimates of the contrasts you have coded. However, that is not the case. Let's look at just one three-level factor. Your coding is $$ \begin{align*} \mu_1 &= \beta_0 -2\beta_1 + 0\beta_2\\ \mu_2 &= \beta_0 +\beta_1 -\beta_2\\ \mu_3 &= \beta_0 + \beta_1 + \beta_2 \end{align*} $$ Note that the contrast codings are weights applied to the coeffiicients, not to the means. You have to solve for $\beta$ to understand what's being estimated:

> X = matrix(c(1,1,1,-2,1,1,0,-1,1), ncol=3)
> X
     [,1] [,2] [,3]
[1,]    1   -2    0
[2,]    1    1   -1
[3,]    1    1    1
> solve(X)
           [,1]       [,2]      [,3]
[1,]  0.3333333  0.3333333 0.3333333
[2,] -0.3333333  0.1666667 0.1666667
[3,]  0.0000000 -0.5000000 0.5000000

So $\beta_0 = (\mu_1 + \mu_2 + \mu_3)/3$, $\beta_1 = (-2\mu_1 + \mu_2 + \mu_3)/6$, and $\beta_3 = (\mu_3 - \mu_2)/2$.

Really, contrast coding is not as straightforward as it may seem. And it gets tougher when you have more than one factor, and perhaps some interactions and covariates. It is far safer to use software that will get it right. If M is the data frame of contrasts shown in the OP,

library(emmeans)
EMM <- emmeans(model, ~ F2 * F1)  # creates the order shown in OP
contrast(EMM, M[, 3:6])

This works correctly whether or not you used those contrast codings in fitting the model.

$\endgroup$
1
  • $\begingroup$ Thank you for the extensive answer. I totally agree that it's by far not as easy as it's made out to be. I think I understood your example from the literature, I've been follwing this tutorial $\endgroup$
    – c06n
    Commented Jun 7, 2021 at 12:00
1
$\begingroup$

Ok after too much head banging I found a couple of issues. Here is hoping this post can save others from the pain.

  1. The estimates are indeed scaled, or rather, not scaled. I followed this tutorial to create my contrasts. It uses the inverse of a hypothesis matrix to generate the contrast matrix. This works well for cases when the compared groups have all the same size. If they differ, like in my example or also when using Helmert contrasts, the estimates are only a fraction of what they should be. Contrasts can / should be divided by the number of groups they calculate estimates for, see here for an overview, and here for the calculation. When using this approach (and not the inverse), the estimates are almost correct ...

  2. Except at first not quite in emmeans::contrasts(). Why? Because the reference group is different from the Helmert contrast-standard, which I expected. It compares the tested group with the mean of all groups involved in the comparison, not just the other groups. It's mean(A + B + C) - C, not mean(A + B) - C. It would be a massive improvement for emmeans to always signify and provide the reference group, because that is absolutely crucial for interpretation of the results.

So why isn't this issue better known? Because the estimates really are just the differences of the compared (reference) groups. The contrast coding 'just' tests for significance, which does not change when scaling the contrasts. Calculating them from the raw data appears to yield the exact same results, so the exact values are of no concern.

I think it's still important to know in order to understand contrast coding.

Btw, not trying to bash emmeans in any way. It's a great package.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.