Group level effects (odds ratio) - dummy coding vs. effect

Question

I ran a logistic model in R's brms, with a categorical variable (condition) with two levels as predictor. When interpreting the results, I realised I don't know how to deal with (the default) dummy coding $(0, 1)$ vs. effect coding $(-1, 1) or (-.5, .5)$ of my categorical variable. I know the intercept is that of condition $0$ for dummy coding, the mean of the two conditions for effect coding. And, as I understand it (I might be wrong here), the effect of the categorical predictor will be equivalent to switching from condition $0$ to $1$ for dummy coding (the intercept + the effect) and will denote the mean deviation from the (mean) intercept for sum coding. For deviation coding, I get the overall difference between the intercept of the two conditions. (At least this is my interpretation of this post.) However, I'm stuck regarding these matters:

Group level effects

How does the coding affect my random intercepts? Do I get the deviation of condition $0$ from its respective intercept with dummy coding? I.e. learn nothing about the variance of condition $1$? And, accordingly, get the overall deviation from the (mean) intercept of both conditions in the case of effect coding?

Odds ratio

Matters are complicated by the fact that my results are in odds ratio. I just tried dummy coding vs. effect coding, but the results slightly differ from what I expected. This will probably get clearer in the example below.

Example and Output

Using brms, I fit a model of the form dropout ~ cond + (1|ID) + (1|marker), where dropout is the dropout rate in each trial (family = bernoulli()) and cond refers to the experimental condition, which could either be "mat" (mathematical) or "nonmat"(non-mathematical). I added group level effects in the form of random intercepts for each participant (ID) and each stimulus (marker). In the column dropout, 1 indicates a trial being considered for analysis, 0 indicates a dropped-out trial.

The model output presented is already converted to odds ratio via exp().

Dummy coding (default)

With the default dummy coding (i.e. brms automatically codes mat = 0 and nonmat = 1), I get the following output:

  .variable            .value .lower .upper .width .point .interval
  <chr>                 <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 b_condnonmat           1.24  0.940   1.67   2.59 mean   hdi      
2 b_Intercept           42.6  24.4    78.1    2.59 mean   hdi      
3 sd_ID__Intercept       4.79  3.00    8.56   2.59 mean   hdi      
4 sd_marker__Intercept   1.38  1.00    1.72   2.59 mean   hdi 

As I understand it, the odds of a trial being analysed (no dropout) are $42.6$ in the "mat" condition (i.e. $42.6 / (1 + 42.6) \approx 97.7 \%$). For the "nonmat" condition, the odds then change by the factor 1.24, meaning that the odds of a trial being analysed are $42.6 * 1.24 = 52.824$ in the "nonmat" condition ($\approx 98.14 \%$). For the group level effects, I'd say that for each participant (ID) the odds of a trial being analysed in the "mat" condition vary, on average, by a factor of $4.79$ - but I'm unsure about that.

Sum coding

To check my assumptions, I coded "mat" as 1 and "nonmat" as -1 and ran the model again, which gives me:

  .variable            .value .lower .upper .width .point .interval
  <chr>                 <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 b_cond                0.898  0.773   1.03   2.59 mean   hdi      
2 b_Intercept          47.7   27.7    86.2    2.59 mean   hdi      
3 sd_ID__Intercept      4.82   2.92    8.23   2.59 mean   hdi      
4 sd_marker__Intercept  1.39   1.01    1.74   2.59 mean   hdi  

Here, my interpretation is that the odds of a trial being analysed, regardless of the experimental condition, are $47.7$ ($\approx 97.95\%$). The condition means a deviation from these odds by the factor of $.898$ and since "mat" was coded as 1, I assumed that I'd get to the odds of a trial being analysed in the "mat" condition with $47.7 * .898 = 42.8346$. However, this is a slight deviation of the $42.6$ I got above. Am I making a mistake here or is this just due to rounding or slight deviations in the model fitting process? (I ran both models with the same seed.) Also, I don't know how to get to the odds for the "nonmat" condition in this case.

The group level effects look slightly different in this model, but not much. I thought that in this model, they refer to the (overall) intercept, so for each participant (ID) the odds of a trial being analysed vary, on average, by a factor of $4.82$.

But since the estimates of the group level effects are so similar, I was wondering if they just happen to be similar for both coding types or if their interpretation is not affected by the coding of the predictor. I was also wondering if the coding changed the model fitting process ever so slightly (even though the seed was the same and it shouldn't matter in principle), which explains the different results.

Deviant coding

Lastly, I coded "mat" as .5 and "nonmat" as -.5 and ran the model again, which gives me:

  .variable            .value .lower .upper .width .point .interval
  <chr>                 <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 b_cond                0.806  0.607   1.06   2.59 mean   hdi      
2 b_Intercept          48.2   27.0    86.2    2.59 mean   hdi      
3 sd_ID__Intercept      4.65   2.89    8.28   2.59 mean   hdi      
4 sd_marker__Intercept  1.40   1.02    1.76   2.59 mean   hdi  

Since the intercept should have the same interpretation for sum and deviant coding, I assume that the mere change of the coding scheme also affects my model fitting process. Are there any guidelines for best practices?

Any help for navigating through log land is much appreciated!

Useful stuff I found

The link mentioned above provides a great overview for the interpretation of different coding schemes - unfortunately, not for random effects. It does, however, mention this source where it gets evident that random effect interpretation does differ for different predictor coding schemes - it's just not clear to me how exactly.

Answer 1

Okay, I think I've figured this one out. First some notation. We have an observation $i$ in group $j$ where the outcome of interest $Y$ follows a Bernoulli distribution and you'll be placing a linear predictor on it's $p$ parameter using the inverse logit link:

\begin{align} Y_{ij} &\sim \text{B}(p_{ij})\\ p_{ij} &= \text{logit}^{-1}(\alpha + \beta c_{ij} + \gamma_j), \end{align}

where $\alpha$ is the global intercept, $\beta$ is the change associated to changes in the condition $c_{ij}$ and $\gamma_j$ is the group-specific intercept.

It can be shown that the specific value we choose for contrasts (e.g. $-\frac{1}{2}$ and $\frac{1}{2}$ or $-1$ and $1$) will only change the scale of $\beta$, so I'll pick the first one because it leaves $\beta$ unchanged compared to dummy coding and hence makes the math easier.

It's easy to check that the choice of coding won't affect the meaning of the slope, which is the difference between each of the conditions in a given group:

\begin{align} (\alpha + \beta + \gamma_j) - (\alpha + \gamma_j) &= \beta\\ (\alpha + 0.5\beta + \gamma_j) - (\alpha - 0.5\beta + \gamma_j) &= \beta \end{align}

The presence of group-specific effects didn't make a difference there. If we had used $1$ instead for the contrasts, the difference for the second equation would be $2\beta$, meaning the model would return a slope that is scaled to half of that under dummy coding.

Now, does this mean there are no differences at all? Not necessarily.

It seems reasonable to ask that, regardless of coding, our model should predict the same values given a group and condition. Let $p_{0j}$ be the predicted probability under dummy coding and $p_{-j}$ the predicted probability under contrast coding. Since we don't yet know what the effect of coding will be on the global and group intercepts, let's also distinguish those: $\alpha_D$ and $\gamma_{Dj}$ will be the parameters under dummy coding while $\alpha_C$ and $\gamma_{Cj}$ the ones for contrast coding. We then have

\begin{align} p_{0j} = p_{-j}\\ \Rightarrow \text{logit}^{-1}(\alpha_D + \gamma_{Dj}) = \text{logit}^{-1}(\alpha_C - 0.5\beta + \gamma_{Cj})\\ \Rightarrow \alpha_D + \gamma_{Dj} = \alpha_C - 0.5\beta + \gamma_{Cj} \end{align}

This relates our dummy coding parameters to the contrast ones. But there's two of them and a single equation, so we need more information to assign them unique values.

For the sake of argument, let's see what happens if we set the group intercepts to be equal:

$$\alpha_D = \alpha_C - 0.5\beta$$

This is the result you're already familiar with, where $\alpha_C$ is bigger than $\alpha_D$ by half a $\beta$, or we can say it's placed at the "average" of both conditions.

If instead we want to keep the global intercepts unchanged, we'd have

$$\gamma_{Dj} = \gamma_{Cj} - 0.5\beta$$

which by itself isn't any big revelation. Until we remember that those group intercepts have a distribution! That distribution is commonly assumed to be a zero-centered normal. But $\beta$ in general won't be zero, so the assumption that $\gamma_{Dj}$ comes from $N(0,\sigma)$ means that $\gamma_{Cj}$ must come from $N(0.5\beta,\sigma)$.

Now, what's going to happen when you try to fit the model? Since a zero-centered normal is the default assumption for every package I've seen, this means the global intercept would be the one to "absorb" the changes introduced by contrast coding.

Now for some simulations

## Simulation function
datasim <- function(npercg, ngroup, alpha, beta, sigma){
  # npercg: number of observations per condition per group
  # ngroup: number of groups
  # beta: effect size for the condition
  # sigma: random intercept spread

  if(npercg<1|ngroup<1|sigma<=0)stop("Something's wrong with your inputs")

  c0 <- rep(0, npercg) 
  c1 <- c0 + 1
  cx <- c(c0, c1)

  gid<- 1:ngroup
  g  <- rnorm(ngroup, mean = 0, sd = sigma)
  gv <- rep(g, each = 2*npercg)

  p  <- 1/(1 + exp(-(alpha + beta*cx + gv)))

  data.frame(y = rbinom(length(p), 1, p), x = cx,
             probability = p, group = as.factor(rep(gid, each = 2*npercg)), rintercept = gv)

}

### Simulate dataset
set.seed(69420)
d1 <- datasim(40, 30, -0.25, 1, 0.33)

### Fitting
library(lme4)

freq_dummy <- glmer(y ~ x + (1|group), data = d1, family = binomial)
freq_contr <- glmer(y ~ I(x-0.5) + (1|group), data = d1, family = binomial)
freq_contr2<- glmer(y ~ I(2*(x-0.5)) + (1|group), data = d1, family = binomial)

library(brms)
options(mc.cores = 4)
cl <- list(max_treedepth = 13)

bayes_dummy <- brm(y ~ x + (1|group), data = d1, family = bernoulli(), control = cl)
bayes_contr <- brm(y ~ I(x-0.5) + (1|group), data = d1, family = bernoulli(), control = cl)
bayes_contr2<- brm(y ~ I(2*(x-0.5)) + (1|group), data = d1, family = bernoulli(), control = cl)

### Results
matrix(c(summary(freq_dummy)$coefficients[,1],
   summary(freq_contr)$coefficients[,1],
       summary(freq_contr2)$coefficients[,1]), nrow = 2)

#            [,1]      [,2]      [,3]
# [1,] -0.2283242 0.2948815 0.2948815
# [2,]  1.0464066 1.0464027 0.5232013

summary(freq_dummy)$varcor
summary(freq_contr)$varcor
summary(freq_contr2)$varcor

# Groups Name        Std.Dev.
# group  (Intercept) 0.33464
# group  (Intercept) 0.33464
# group  (Intercept) 0.33464

matrix(c(summary(bayes_dummy)$fixed[,1],
     summary(bayes_contr)$fixed[,1],
         summary(bayes_contr2)$fixed[,1]), nrow = 2)

#            [,1]      [,2]      [,3]
# [1,] -0.2312706 0.2963076 0.2960156
# [2,]  1.0481774 1.0492645 0.5235235

summary(bayes_dummy)$random
summary(bayes_contr)$random
summary(bayes_contr2)$random

#                Estimate  Est.Error  l-95% CI u-95% CI Eff.Sample     Rhat
# sd(Intercept) 0.3616682 0.07577461 0.2328071 0.524049   1281.155 1.001547
# sd(Intercept)   0.36224 0.07594461 0.2316291 0.5347441  1396.271 1.001551
# sd(Intercept) 0.3632387 0.07466388 0.2350172 0.5263185  1571.552 1.001344

I'd say this confirms the reasoning we developed previously. As for the small differences in the Bayesian models, how can we know they are are solely attributable to Montcaro error? It's pretty simple: just run the same model a couple of times and you'll notice that the same-model differences are enough to make the inter-model differences overlap zero. I didn't show the results for that, but I did it a few times and feel pretty confident that the differences are just MC error.