Interpretation of regression coefficients with multiple categorical predictors

Question

This question has an UPDATE.

There is a nice answer HERE regarding how to interpret regression coefficients when predictors each consist of two categories in R. But imagine we have students' sex (boys, girls) and the school-gender system (boy-only, girl-only, mixed) in a model like: y ~ sex + schoolgend.

Here, we expect 4 coefficients. I wonder how to interpret the third coefficient (+.175 see below)?

Say, the (Intercept) represents the boys mean in mixed schools:

               Estimate
(Intercept)      -0.189
schgendboy-only   0.180
schgendgirl-only  0.175 <-- This one!
sexgirls          0.168

My interpretations of the coefficients are as follows:

     (intercept): mean of y for boys in mixed schools = -.189

 schgendboy-only: diff. bet. boys in boy-only vs. mixed schools = +.180

schgendgirl-only: diff. bet. ?????????????????????????????????? = +.175

        sexgirls: diff. bet. girls vs. boys in mixed schools = +.168

If my interpretation logic for all other coefs is correct, then, this third coef. (schgendgirl-only) must mean:

diff. bet. boys in girl-only vs. mixed schools = +.175! (which makes no sense!)

ps. I know I will end-up interpreting +1.75 as: diff. bet. girls in girl-only vs. mixed schools BUT this doesn't follow the interpretation logic for other coefs PLUS there are no labels in the output to show what's what!

UPDATE

It seems that schgendgirl-only can meaninglessly represent: diff. bet. boys in girl-only vs. mixed schools just like it can meaningfully represent: diff. bet. girls in girl-only vs. mixed schools

The reason is that the two coefficients are exactly equal. Here is a reproducible R demonstration (results are shown to 4 decimal places):

library(R2MLwiN) # For the dataset
library(lme4)    # For Model fitting
library(emmeans) # For pairwise contrasts

data("tutorial")

Form <- normexam ~ 1 + standlrt + schgend + sex + (standlrt | school)
model <- lmer(Form, data = tutorial, REML = FALSE)

emmeans(model, pairwise~schgend+sex)$contrast

#contrast                     estimate     SE  df z.ratio p.value
#mixedsch boy - boysch boy    -0.17986 0.0991 Inf -1.814  0.4565
#mixedsch boy - girlsch boy   -0.17482 0.0788 Inf -2.219  0.2287<--This coef. equals
#mixedsch boy - mixedsch girl -0.16826 0.0338 Inf -4.975  <.0001
#mixedsch boy - boysch girl   -0.34813 0.1096 Inf -3.178  0.0186
#mixedsch boy - girlsch girl  -0.34308 0.0780 Inf -4.396  0.0002
#boysch boy - girlsch boy      0.00505 0.1110 Inf  0.045  1.0000
#boysch boy - mixedsch girl    0.01160 0.0997 Inf  0.116  1.0000
#boysch boy - boysch girl     -0.16826 0.0338 Inf -4.975  <.0001
#boysch boy - girlsch girl    -0.16322 0.1058 Inf -1.543  0.6361
#girlsch boy - mixedsch girl   0.00656 0.0928 Inf  0.071  1.0000
#girlsch boy - boysch girl    -0.17331 0.1255 Inf -1.381  0.7388
#girlsch boy - girlsch girl   -0.16826 0.0338 Inf -4.975  <.0001
#mixedsch girl - boysch girl  -0.17986 0.0991 Inf -1.814  0.4565
#mixedsch girl - girlsch girl -0.17482 0.0788 Inf -2.219  0.2287<--This coef.
#boysch girl - girlsch girl    0.00505 0.1110 Inf  0.045  1.0000

Answer 1

I think your question fundamentally goes back to the difference between the interpretation of main-effects in main-effect-only (y ~ cat1 + cat2) vs. interactive (y ~ cat1 * cat2) regression models.

When your model is of the form: y ~ sex + schoolgend. Then, YES, .175 has no choice but to represent both the:

diff. bet. girls in girl-only vs. mixed schools

as well as:

diff. bet. boys in girl-only vs. mixed schools (we'll see why this substantively doesn't make sense in your case)

In main-effect-only models, this interpretation arises from the fact that you can keep any one predictor constant at any of its levels (doesn't matter which one) while varying the other predictor's levels.

So, for schgendgirl-only coef., you're saying holding sex constant at any of its levels (either boys, OR girls) which equivalently means regardless of students' sex; boys or girls, how much y will change if we change schgend from mixed (its reference level) to schgendgirl-only (the level of the tabled coef.)?

In interactive models, however, the main effect of one predictor is measurable ONLY when we keep the other predictors at their reference values/levels.

So, your problem was that you applied the logic of interpretation for the main effects in interactive models to that in main-effect-only models.

Indeed, it makes good sense in your case to let students' sex and schoolgend to interact. Because, as researchers, we may reasonably want to know how y changes for different genders situated in same vs. different school gender systems (a contextual effect).

Now let's see why some of the comparisons based on your model are substantively meaningless. By the way your two categorical predictors are set-up, you're telling your model that sex has two levels, and schgend has three levels. But in reality because boys can't be in girl-only, and girls can't be in boys-only schools (data for such combinations don't exist), you get several irrelevant extrapolations as shown in your post-hoc analysis for boys in girl-only and girls in boy-only schools as you don't let sex and schgend to interact.

So, what to do?

First, recode your schgend to have its real two levels (same_sex [all girls or all boys] vs. mixed schools):

tutorial2<-transform(tutorial,schgend2=ifelse(schgend!="mixedsch","same_sex","mixedsch"))

Second, change your model to be interactive:

Form2 <- normexam ~ 1 + standlrt + schgend2 * sex + (standlrt | school)
model2 <- lmer(Form2, data = tutorial2, REML = FALSE)
round(coef(summary(model2)),4)

                         Estimate Std. Error t value
(Intercept)               -0.1888     0.0514 -3.6767
standlrt                   0.5544     0.0199 27.8071
schgend2same_sex           0.1799     0.0991  1.8141
sexgirl                    0.1683     0.0338  4.9750
schgend2same_sex:sexgirl  -0.0050     0.1110 -0.0454

Notice that the meaning of schgend2same_sex ($0.1799$) now matches your interactive-model interpretation logic i.e., it ONLY represents: diff. bet. boys in boy-only vs. mixed schools.

The corresponding difference for girls is a tiny bit (~$-.005$) smaller. That is, diff. bet. girls in girl-only vs. mixed schools is ~$0.1749$ (you'll see this in the post-hoc table below shown as $.1748$ due to rounding).

Third, perform your desired post-hocs:

emmeans(model, pairwise~schgend2*sex)$contrast

# contrast                      estimate     SE  df z.ratio p.value
# mixedsch boy - same_sex boy    -0.1799 0.0991 Inf -1.814  0.2666 
# mixedsch boy - mixedsch girl   -0.1683 0.0338 Inf -4.975  <.0001 
# mixedsch boy - same_sex girl   -0.3431 0.0780 Inf -4.396  0.0001 
# same_sex boy - mixedsch girl    0.0116 0.0997 Inf  0.116  0.9994 
# same_sex boy - same_sex girl   -0.1632 0.1058 Inf -1.543  0.4115 
# mixedsch girl - same_sex girl  -0.1748 0.0788 Inf -2.219  0.1180

Now, you get all six possible caparisons between your categorical predictors.

Answer 2

It's the difference between the (predicted) mean of girls in girls-only schools and the (predicted) mean of girls in mixed schools. You can see this by looking at the design matrix or solving for a predicted value using the fitted regression equation. Let's make some simple data and work through this.

d = expand.grid(sex=c("boys","girls"), 
                schoolgend=c("boy-only", "girl-only", "mixed"))
d   = d[-c(2,3),]
d$y = c(-0.189+0.180, -0.189+0.168+0.175, -0.189, -0.189+0.168)
d   = as.data.frame(lapply(d, rep, 2))
d   = d[with(d, order(schoolgend, sex)),]
d$y = d$y + rep(c(-.001, .001), times=4)
d$schoolgend = relevel(d$schoolgend, ref="mixed");    d
aggregate(y~sex+schoolgend, d, mean)
#     sex schoolgend            y
# 1  boys      mixed -0.189
# 2 girls      mixed -0.021
# 3  boys   boy-only -0.009
# 4 girls  girl-only  0.154
model.matrix(lm(y~sex+schoolgend, d))  # this is the design matrix
#   (Intercept) sexgirls schoolgendboy-only schoolgendgirl-only
# 1           1        0                  1                   0
# 5           1        0                  1                   0
# 2           1        1                  0                   1
# 6           1        1                  0                   1
# 3           1        0                  0                   0
# 7           1        0                  0                   0
# 4           1        1                  0                   0
# 8           1        1                  0                   0
m = lm(y~sex+schoolgend, d)
t(t(coef(m)))
# (Intercept)         -0.189
# sexgirls             0.168
# schoolgendboy-only   0.180
# schoolgendgirl-only  0.175

Looking at the design matrix, it's easy to see that schoolgendgirls-only marks off the difference between girls in mixed schools and girls in girls-only schools. You can also see this by looking at the complete regression equation.

\begin{align} \hat{y} &= -0.189 + 0.168 x{\rm(sex=girls)} + 0.180 x{\rm(sg=boyonly)} + 0.175x{\rm (sg=girlonly)} \\ \hat{y} &= -0.189 + 0.168 \cdot 1 + 0.180 \cdot 0 + 0.175 \cdot 1 \\ \hat{y} &= -0.189 + 0.168 + 0.175 \end{align}

---

Update:

The main question appears to be whether the interpretation below is correct.

I have learned to interpret any main effect coef for a categorical predictor by thinking of that coef. as something that can differ from its reference category to affect "y" holding any other categorical predictor in the model at its reference category

This is not quite correct, although it is a common (and intuitive) way to explain it (and which I have used myself). The truth is that is the difference with its reference level holding all other variables constant. The other variables can be held constant at their means / reference levels, or at any other value. It may help to read my answer here: What does "all else equal" mean in multiple regression?

Answer 3

The answers so far address interpretations that implicitly ignore the multilevel nature of the model. They are all correct though, that you don't have interactions, so the interpretations hinge on the tricky issue of boys-only or girls-only schools not having the other sex. The R code is helpful but let's put this in terms of the coefficients and see where we're at:

normexam ~ 1 + standlrt + schgend + sex + (standlrt | school)

Level 1 (student):

$$exam_{ij} = \beta_{0j} + \beta_{1j}standlrt_{ij} + \beta_{2j}sex_{ij} + r_{ij}$$

Level 2 (school):

$$\beta_{0j} = \gamma_{00} + \gamma_{01}schgendboys_{j} + \gamma_{02}schgendgirls_{j} + u_{0j}$$ $$\beta_{1j} = \gamma_{10} + u_{1j}$$ $$\beta_{2j} = \gamma_{20}$$

Combined: $$exam_{ij} = \gamma_{00} + \gamma_{01}schgendboys_{j} + \gamma_{02}schgendgirls{j} + \gamma_{10}standlrt{ij} + \gamma_{20}sex_{ij} + u_{0j} + u_{1j}standlrt_{ij} + r_{ij}$$

Let's hold sex constant at male. In this model (with dummy coding), $\gamma_{00}$ (the intercept) is the grand mean exam score in mixed schools. $\gamma_{01}$ is the difference in school average between mixed and boys-only, $\gamma_{02}$, difference in school average between mixed and girls only, $\gamma_{10}$ is the pooled estimate (average) standlrt slope across all schools, $\gamma_{20}$ is the overall difference between male and female students (non-random).

It's easier to think about it in the following way: within each school we're supposedly estimating the standlrt slopes controlling for the sex of the student, and then between schools, we're examining the effect of school type (mixed, boys-only, girls-only) on the school's average exam scores, while controlling for the level-1 variables within a school. Within the single sex schools, obviously the effect of sex cannot be estimated, but we still get a coefficient $\gamma_{20}$ because we're actually estimating it using the entire sample, and not letting it vary across schools. This works because in arriving at MLM estimates, the individual school regressions are actually not performed. Instead, estimates of variance components (the r's and the u's) enable the fitting of the combined model. One can then obtain the school-level intercepts and slopes (you should try that to see if some are NA's).

Hope this helps!