2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Small pet peeve: I believe you are referring to the same coding (-0.5, 0.5) by sum, effect, deviant and deviation coding. Please try to keep the terminology consistent! I'm interested in an answer to this myself, I'll do some simulations and report back. $\endgroup$ – zipzapboing Aug 22 at 2:05
  • $\begingroup$ Maybe I wasn't clear in my formulation (will check that), but based on this link I thought the following to be the case: effect coding = sum or deviation coding. While sum coding (-1, 1) and deviation coding (-.5, .5). $\endgroup$ – einGlasRotwein Aug 22 at 9:04
  • 1
    $\begingroup$ Apart from the difference between dummy coding vs effect coding, you should be aware that the interpretation of the fixed effects coefficients is conditional on the random effects. For more on this, check here: stats.stackexchange.com/questions/365907/… $\endgroup$ – Dimitris Rizopoulos Aug 24 at 6:06
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks a lot! For the math and especially the simulations. Just to confirm that I understood that correctly: The coding scheme doesn't change the estimate of the group-level effects, but when I use contrast coding, it messes up the global intercept because the center of the distribution of group level intercepts isn't at 0 anymore? That means, whenever I use group level effects with contrast coding, my global intercepts will be off? That would be quite important because a lot of papers I've seen use contrast coding because it usually makes the global effects more "intuitive". $\endgroup$ – einGlasRotwein Sep 4 at 17:47
  • $\begingroup$ Hmm, wait, you mention group-level intercepts and a global intercept, which I showed above, but then also mention group-level effects and global effects. Can you clarify what you mean with these terms? In my example, and in mixed models in general, the "slope" ($\beta$ here) is interpreted as the effect conditional on group membership, which one could parse as "group level effect", however it's a single identical parameter shared across all groups, so you could also say it's a "global effect". $\endgroup$ – zipzapboing Sep 5 at 3:27
  • $\begingroup$ Anyway, I definitely should have been more concise here, sorry about that! I was discovering the answer as I wrote it. To summarize: contrast coding does not "mess up" the global intercept in the sense of making it take some undesirable value, but it does change its value and and its interpretation. My point about the group level intercepts is that they could be what changes instead of the global intercept, but because all statistical software forces their mean to be zero, then you're also guaranteed that will not happen and the global intercept will always be the only thing which changes. $\endgroup$ – zipzapboing Sep 5 at 3:34
  • $\begingroup$ One small pedantic point: that the model doesn't change with coding isn't always true; if you use a strong prior on your global intercept, then you would be restricting how much it can change and so a diference in coding would actually lead to somewhat different predictions. But this is why Bayesian modeling guidelines strongly recommend that global intercept terms are given wide priors and brms follows that advice so you shoud not run into any issues related to this in practice. $\endgroup$ – zipzapboing Sep 5 at 3:43
  • $\begingroup$ Sorry, minefield of terms that one can get wrong. I learned "fixed effects" as population level effects and "random effects" as group level effects, i.e. there are different sub groups in the data (other than e.g. the experimental condition, which I'd never refer as "group effect" to). What I'm now confused about is: The group intercept tells me how much the global intercept varies per subgroup in my data. With dummy coding, that intercept will be a different one than with contrast coding, and I thought this is what you were referring to. $\endgroup$ – einGlasRotwein Sep 6 at 4:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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