1

EDIT: I think I have mistaken the names of the coding systems, so I changed it (in bold). The content has not changed at all, though, so I would still appreciate any answer. END EDIT

I'm running linear generalized mixed effects regression (glmer) on binomial data. To simplify: I have two categorical predictors, A and B. A has two levels (a1 and a2) while B has three (b1, b2, b3). I used a helmert coding system for B, i.e. one contrast codes for the difference between b2 and b1, and the second contrast codes for the difference between b3 and the average of b1 and b2. A is sum coded. With these coding systems, I get that the interaction of A and the first contrast of B (b2>b1) is significant at p=0.08 (see results below). When I change the coding system of B to difference coding (i.e. adjacent levels are compared; contrast I: b2>b1, contrast II: b3>b2), I get that the interaction of A with the first contrast (b2>b1) is significant at p=0.04 (see results below). How can this be? The first contrast of B is the same in both coding systems.

The contrast matrix I used for helmert coding is the following:

> helmert.codes
    b2 vs b1   b3 vs mean(b2,b1)
b1     -0.5     -0.3333333
b2      0.5     -0.3333333
b3      0.0      0.6666667

And the contrast matrix for difference coding was this:

> difference.codes
      b2 vs b1  b3 vs b2
b1       -1        0
b2        1       -1
b3        0        1

After searching the forum, I ran across this post in which @barnhillec explains that for difference coding, it is the pseudo-inverse of the above matrix that should be used. Based on this, I tried to use the following matrix instead:

> inverse.difference.codes
     b2 vs b1     b3 vs b2
b1 -0.6666667   -0.3333333
b2  0.3333333   -0.3333333
b3  0.3333333    0.6666667

With this matrix, I once again get that the interaction of A and the first contrast of B is significant at the level of p=0.08 (although there are still some differences three places after the decimal).

So my questions are:

  1. Which matrix should I use if I want a difference coding system? Assuming the inverse is the one, could someone please help me put sense into it? I don't understand how the inverse is the one that codes for these contrasts.
  2. Should I expect any difference between the two coding systems (helmert vs difference) in the first contrast of B (and its interactions), at all?

* LATER ADDITION: RESULTS *

Results with helmert coding (relevant effects only):

Random effects:
 Groups      Name                 Variance Std.Dev. Corr                   
 Participant (Intercept)          0.50777  0.7126                          
             A                    0.22049  0.4696   -0.19                  
             Bb2_vs_b1            0.10382  0.3222   -0.82  0.70            
             Bb3_vs_mean(b2,b1)   0.10101  0.3178   -0.75  0.17  0.59      
             A:Bb2_vs_b1          0.31710  0.5631   -0.26 -0.07  0.11  0.82
             A:Bb3_vs_mean(b2,b1) 0.04621  0.2150    0.08  0.07 -0.14  0.01  -0.08

Fixed effects:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)            3.22267    0.17793  18.112  < 2e-16 ***
A                     -0.76855    0.13976  -5.499 3.82e-08 ***
Bb2_vs_b1             -1.47336    0.27862  -5.288 1.24e-07 ***
Bb3_vs_mean(b2,b1)    -1.31688    0.19104  -6.893 5.46e-12 ***
A:Bb2_vs_b1           -0.52065    0.29891  -1.742 0.081545 .  
A:Bb3_vs_mean(b2,b1)  -0.29022    0.18209  -1.594 0.110987

Results with difference coding, matrix difference.codes:

Random effects:
 Groups      Name           Variance Std.Dev. Corr                         
 Participant (Intercept)    0.50777  0.7126                                
             A              0.22049  0.4696   -0.19                        
             Bb2_vs_b1      0.22922  0.4788   -0.88  0.55                  
             Bb3_vs_b2      0.17957  0.4238   -0.75  0.17  0.84            
             A:Bb2_vs_b1    0.32473  0.5699   -0.24 -0.05  0.41  0.81      
             A:Bb3_vs_b2    0.08216  0.2866    0.08  0.07 -0.09  0.01  0.17

Fixed effects:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)         3.22267    0.17773  18.133  < 2e-16 ***
A                  -0.76855    0.13975  -5.499 3.81e-08 ***
Bb2_vs_b1          -2.35130    0.34544  -6.807 9.99e-12 ***
Bb3_vs_b2          -1.75584    0.25423  -6.907 4.96e-12 ***
A:Bb2_vs_b1        -0.71409    0.35412  -2.017 0.043746 *  
A:Bb3_vs_b2        -0.38696    0.24241  -1.596 0.110415

Results with difference coding, matrix inverse.difference.codes:

Random effects:
 Groups      Name           Variance Std.Dev. Corr                         
 Participant (Intercept)    0.50777  0.7126                                
             A              0.22049  0.4696   -0.19                        
             Bb2_vs_b1      0.10383  0.3222   -0.82  0.70                  
             Bb3_vs_b2      0.06657  0.2580   -0.41 -0.23  0.10            
             A:Bb2_vs_b1    0.31709  0.5631   -0.26 -0.07  0.11  0.94      
             A:Bb3_vs_b2    0.13516  0.3676    0.24  0.09 -0.17 -0.66 -0.81

Fixed effects:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)         3.22267    0.17779  18.126  < 2e-16 ***
A                  -0.76856    0.13973  -5.500 3.79e-08 ***
Bb2_vs_b1          -1.47342    0.27861  -5.288 1.23e-07 ***
Bb3_vs_b2          -0.58017    0.19069  -3.042 0.002347 ** 
A:Bb2_vs_b1        -0.52054    0.29897  -1.741 0.081661 .  
A:Bb3_vs_b2        -0.02994    0.19775  -0.151 0.879669
Galit
  • 177
  • 1
  • 10
  • 2
    for what it's worth, this isn't (AFAICT) specific to mixed models at all, it's most likely generic to linear model contrasts. Can you show us the coefficient tables for both versions? That might help to explain the results. – Ben Bolker Dec 31 '18 at 19:44
  • I added now all three models' estimates. – Galit Dec 31 '18 at 20:32
  • I think instead of worrying about coding systems, you should tell us what it is you are trying to accomplish in your analysis. Are you trying to compare? For example, are you trying to compare levels of a variable with the mean of the subsequent levels of the variable or perhaps you want to Compares each level to the reference level or comparing each level to an overall grand mean? Depending on your answer, we can guide you to which coding scheme will work best for you. – StatsStudent Jan 02 '19 at 20:38
  • Basically, I want to compare adjacent levels. I was using at first helmert coding because it is orthogonal, but then I switched to difference coding because it better suits the question. I did not expect it to change the estimate of the first contrast, because it is the same contrast (b2 vs b1). Thus, there must be something here that I did wrong in defining the matrices, which I don't understand. – Galit Jan 02 '19 at 21:37

1 Answers1

2

After some web research, I finally found an answer. So I am answering my own question and hopefuly this will be helpful for future similar issues. I base my answer on this great explanation by Bill Venables. It goes like this:

There is a difference between a coding matrix and a contrast matrix. A coding matrix defines the means in terms of the betas, while a contrast matrix defines the betas in terms of the means. Therefore, one is simply the inverse of the other. The matrices in my question (helmert.codes and inverse.difference.codes) are coding matrices (though they are the input of the contrasts() function in R, quite confusing). This means that they define means based on betas, and not the other way around. Names of rows and columns in my example were misleading. Since this is a transformation matrix (from betas to means), the columns are the betas. For example, take the matrix inverse.difference.codes, re-written here (with a b0 column):

> inverse.difference.codes # (coding matrix)
    b0        b1           b2
m1   1    -0.6666667   -0.3333333
m2   1     0.3333333   -0.3333333
m3   1     0.3333333    0.6666667

Each row defines the linear combination of betas (i.e. the explicit equations of the model): \begin{equation} \mu_1=\beta_0-\frac{2}{3}\beta_1-\frac{1}{3}\beta_2 \end{equation} \begin{equation} \mu_2=\beta_0+\frac{1}{3}\beta_1-\frac{1}{3}\beta_2 \end{equation} \begin{equation} \mu_3=\beta_0+\frac{1}{3}\beta_1+\frac{2}{3}\beta_2 \end{equation} If we want to solve this set of equations and find what each beta represents (i.e. taking the inverse of the coding matrix), we get the following solution: \begin{equation} \beta_0 = \frac{1}{3}\mu_1+\frac{1}{3}\mu_2+\frac{1}{3}\mu_3 \end{equation} \begin{equation} \beta_1=\mu_2-\mu_1 \end{equation} \begin{equation} \beta_2=\mu_3-\mu_2 \end{equation} And in matrix notation, this is the transpose of the contrast matrix I called difference.codes, re-written here (with the intercept b0, and with accurate row and column names):

> difference.codes # (contrast matrix)
            b0        b1       b2
m1     0.3333333      -1        0
m2     0.3333333       1       -1
m3     0.3333333       0        1

To summarize, the contrasts() function gets as input a coding matrix, without an intercept (no b0) column, and not a contrast matrix. Therefore, if I want my betas to code for differences between adjacent levels, I should indeed use the following matrix:

> inverse.difference.codes # (coding matrix)
     b1 (m2-m1)   b2 (m3-m2)
m1  -0.6666667   -0.3333333
m2   0.3333333   -0.3333333
m3   0.3333333    0.6666667

First beta is the contrast between m2 and m1, second beta is the contrast between m3 and m2.

Galit
  • 177
  • 1
  • 10