6
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ Please clarify what you mean by the "(global) population mean." Are the four levels of the categorical variable equally represented both in your data sample and in the population? Otherwise, a "global population mean" could be a good deal different from, say, the mean of the means of the 4 groups. $\endgroup$ – EdM Sep 26 at 15:55
  • $\begingroup$ @EdM By global mean I mean the expectation of $Y$. The categorical variable $\tilde{x}$ may be assumed to be multinomially distributed with unequal class probabilities. $\endgroup$ – tomka Sep 26 at 15:58
  • $\begingroup$ Specifically, the expectation of $Y$ in your data sample, perhaps with the assumption that the group representation in your sample represents that in the population? $\endgroup$ – EdM Sep 26 at 16:00
  • $\begingroup$ Please look through this answer with examples stats.stackexchange.com/a/221868/3277. You probably want deviation (effect) type of coding of your factor. $\endgroup$ – ttnphns Sep 26 at 16:02
  • $\begingroup$ @EdM The estimated intercept of the model should be equal to the sample mean, and the sample mean estimates the population mean. See my edit for an example with two categories. $\endgroup$ – tomka Sep 26 at 16:07
9
$\begingroup$

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

enter image description here


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 
| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ also I think referred to as "sum-to-zero" contrasts, or "effect coding" stats.idre.ucla.edu/other/mult-pkg/faq/general/… ... For what it's worth, your answer only works to retrieve the grand mean if the original data set is balanced ... $\endgroup$ – Ben Bolker Sep 26 at 23:46
  • 2
    $\begingroup$ @BenBolker Agreed. If unbalanced it wil be a weighted average of the means in each group, I believe ? $\endgroup$ – Robert Long Sep 27 at 6:05
  • 1
    $\begingroup$ Yes, that's right. In principle there should be a way to "counter-weight" the coded values (probably by just multiplying or dividing each indicator by the population fraction), but this adjustment is usually done post hoc (computing least-squares or expected marginal means) $\endgroup$ – Ben Bolker Sep 27 at 12:17
  • 4
    $\begingroup$ The improvement process is generally one in response to comments on the original question, not on comments to an answer that you don't like. Anyway I won't comment further. I tried my best to answer the question posed and have offered my advise on how you might be able to solve the problem of unbalanced data which was not apparent from the question originally asked. I would be happy to try to answer a new question. $\endgroup$ – Robert Long Sep 28 at 10:15
  • 2
    $\begingroup$ OK fair enough, I've updated my answer :) $\endgroup$ – Robert Long Sep 28 at 10:59

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.