Which categorical variable encoding should I use if I want the model intercept to have the interpretation of the global mean?

Question

In a mixed effects model

$$ y_{ij} = \beta_{00} + \beta_{01} x_{1i} + \beta_{02} x_{2i} + \beta_{03} x_{3i} + u_i + \epsilon_{ij}$$

where $x_1, x_2, x_3$ are dummy variables coding the levels of a discrete (multinomial) variable $\tilde{x}$ that has more than two levels (here four), I want to give the intercept $\beta_{00}$ the interpretation of the (global) population mean, which is $E(y_{ij})$.

Right now the covariates are coded in a way that the intercept is interpreted as the mean of the reference category of $\tilde{x}$.

Is there a way to achieve this?

I have found a good overview on effect coding, but this type of mean coding is not part of it.

Edit: I just remembered how to do this for a variable $\tilde{x}$ that has two categories only. Then we have the model

$$ y_{ij} = \beta_{00} + \beta_{01} x_{1i} + u_i + \epsilon_{ij}$$

where the Dummy $ x_{1i}$ is defined to be $(1-p)$ if $\tilde{x} =1$ and it is $(-p)$ if $\tilde{x} =0$, where $p$ is the proportion with $\tilde{x}=1$.

Edit 2: Following the reply by Robert Long, deviation coding can be used when the number of observations for each level of $\tilde{x}$ are the same. However I am looking for a solution for multinomial $\tilde{x}$ possible with unequal class probabilities. Here is some code to implement deviation coding with multinomial $\tilde{x}$ demonstrating that this coding does not estimate the global mean. I suspect some sort of category weighting on the deviation coded dummies is needed instead (like I did for the two-category case above).

# Code to assess deviation coding for multinomial $xt$    
library(MASS)
library(dplyr)
n = 1000
set.seed(13)
xt = rmultinom(n, 1, c(1/3, 1/3, 1/3))
xt = as.factor( apply( t(t(xt) * c(1,2,3)), 2, sum) )
X <- model.matrix(~ xt)
betas <- c(3, 1, 2)
Y <- X %*% betas + rnorm(n)
mean(Y)

lm(Y ~ xt) %>% coef()   # default treatment coding

contrasts(xt) <- contr.sum(3) # specify deviation coding
lm(Y ~ xt) %>% coef()

Edit 3: Originally the question was titled "Which effect coding (categorical encoding) to use if I want the model intercept to have the interpretation of the global mean?" The title wrongfully suggested that my objective could be reached by effect coding alone. The answer by Robert Long applies to balanced categories and then deviation coding should be used.

Answer 1

If the data are balanced, then deviation coding should work.

Let's look at a simple example:

set.seed(13)
dt <- expand.grid(X1 = LETTERS[1:3], reps = 1:5)
X <- model.matrix(~ X1, dt)
betas <- c(3, 1, 2)
dt$Y <- X %*% betas + rnorm(nrow(dt))
mean(dt$Y)

[1] 4.11413

So we would like the intercept to be 4.11413

If we fit the model with default coding we get:

lm(Y ~ X1, dt) %>% coef()   # default treatment coding

(Intercept)         X1B         X1C 
  3.3430627   0.2867999   2.0264018 

But now if we use deviating coding we get

contrasts(dt$X1) <- contr.sum(3) # specify deviation coding
lm(Y ~ X1, dt) %>% coef()

(Intercept)         X11         X12 
  4.1141299  -0.7710672  -0.4842673 

If the data are unbalanced then you will need to do some post-hoc adjustement.

---

Edit: To address what to do when the data are unbalanced.

In this case, it is easier to work with default treatment coding rather than deviation coding:

> set.seed(1)
> dt1 <- expand.grid(X1 = LETTERS[1:1], reps = 1:5)
> dt2 <- expand.grid(X1 = LETTERS[2:2], reps = 1:3)
> dt3 <- expand.grid(X1 = LETTERS[3:3], reps = 1:2)
> dt <- rbind(dt1, dt2, dt3)
> table(dt$X1)

A B C 
5 3 2 

So the groups are unbalanced.

> X <- model.matrix(~ X1, dt)
> betas <- c(2, 3, 1)
> dt$Y <- 4 + X %*% betas + rnorm(nrow(dt), 0, 1)
> mean(dt$Y)

[1] 7.232203

So we would like to recove 7.23 with a post hoc calculation, which can be acheived fairly easily with

> coef(lm(Y ~ X1, dt))[1] + betas[2] * table(dt$X1)[2]/nrow(dt) + betas[3] * table(dt$X1)[3]/nrow(dt)

(Intercept) 
   7.22927 

Note that the result is not exact due to the combination of imbalance in the groups and the random error. As the error approaches zero, the result becomes exact. Even with error, the result is also unbiased, as we can see from a monte carlo simulation:

n.sim <- 1000
vec.sim <- numeric(n.sim)

for (i in 1:n.sim) {
  
  set.seed(i)

  dt$Y <- 4 + X %*% betas + rnorm(nrow(dt), 0, 1)

  vec.sim[i] <- mean(dt$Y) - (coef(lm(Y ~ X1, dt))[1] + betas[2] * table(dt$X1)[2]/nrow(dt) + betas[3] * table(dt$X1)[3]/nrow(dt))

}

hist(vec.sim)
mean(vec.sim)

[1] -0.003418483

---

Edit: As noted in the comments, we should really be using the coefficient estimates from the model, and doing so will then make the calculation exact:

> coef(lm(Y ~ X1, dt))[1] + coef(lm(Y ~ X1, dt))[2] * table(dt$X1)[2]/nrow(dt) + coef(lm(Y ~ X1, dt))[3] * table(dt$X1)[3]/nrow(dt)
(Intercept) 
   7.232203