Lme4 Crossed Effects Mixed Model w/ Interaction Terms

Question

Greetings Statisticians!

I am trying to create a crossed effects mixed model w/ interaction terms in lme4 to describe a regional analysis of MRI data in mouse brains w/ genotype, treatment, and time as categorical fixed effects to predict MRI signals as the dependent variable. I wish to model these fixed effects as 2x sets of interactions:

1. Treatment x Region
2. Genotype x Region

For these interactions, I have a base category for each variable which I expect to have excluded from the coefficients in the model results - I have successfully achieved this result in earlier models for this project that only utilized main effects in the absence of any interactions. In the case of the brain Region variable, this category is the calculated mean of ALL regions. Random effects include correlated intercepts & slopes for brain region and time, each with mouse subject.

The final model specification is as follows:

formula <- "Y ~ Region:Tx + Region:Transgenic + Week + (Week | Mouse) + 
                (Region | Mouse)"

Since I have 13 regions (+1 "total brain" region), 4x treatment groups, and 2x genotypes, I would expect to have a total of 54x fixed effects coefficients:

13x3 region:treatment coeffs + 13x1 region:genotype coeffs + 1x time coeff + 1x intercept

However, the resulting model is including coefficents that contain the base categories I indicated for a total of 71 coefficients, which is less than the total of 86x fixed effects coefficients I would expect to see if ALL categories were included:

14x4 region:treatment coeffs + 14x2 region:genotype coeffs + 1x time coeff + 1x intercept

As a side note, I'm receiving the following warning message as well:

fixed-effect model matrix is rank deficient so dropping 1 column /   
coefficient
boundary (singular) fit: see ?isSingular

I specify my base categories with the following code:

data <- within(data, Region <- relevel(Region,ref = "Brain"))
data <- within(data, Tx <- relevel(Tx,ref = "Saline"))
data <- within(data, Transgenic <- relevel(Transgenic,ref = "WT"))

Why are my base categories not excluded from the interaction model as they typically are in main effects models?

Answer 1

The statement: "fixed-effect model matrix is rank deficient so dropping 1 column / coefficient" is not a side note. It means that you haven't properly set up your model to work with your data. Some one of your fixed effects is a linear combination of your other fixed effects. You say:

In the case of the brain Region variable, this [reference] category is the calculated mean of ALL regions.

That necessarily makes the values of reference Region category linearly dependent on the values all the individual regions, which then carries through to the interaction terms involving it.

Also, you didn't specify a fixed effect for Region in your formula. I'm not sure how lme4 deals with a random slope that isn't also specified as a fixed effect.

If you want your coefficients for Region to represent differences from a mean value, you could use deviation contrasts for that predictor. Such contrasts are set up by contr.sum() in R.

Then the intercept is the mean among the Region means (not the same as the grand mean if numbers of observations differ), and you get 12 further linearly independent coefficients for 13 regions. The 13th coefficient wouldn't be reported directly by the model, as it is linearly dependent on the Intercept and the other 12 coefficients. The sum of all deviations from the mean must be 0, however, so the 13th deviation is simply the negative of the sum of the other 12. Isabella Ghement illustrates how to work with this coding in detail in this answer. You can get standard errors for the 13th Region by using the coefficient covariance matrix and the formula for the variance of a weighted sum of correlated variables.

Once you remove the linear dependence among your predictors, the model will be fundamentally the same regardless of your predictor-coding scheme. Rather than play with predictor-coding schemes, it might be simpler to use software tools that can provide point estimates and standard errors from a model for any combination of predictor values. The emmeans package is one useful choice for mixed models.

In response to comments (and further thought):

The way that you specified interaction terms without main effects is one source of your problem. Here's a simple example, with 2 categorical predictors, one having 2 levels and another having 3.

A <- factor(rep(letters[1:2],3))
B <- factor(rep(letters[1:3],each=2))

Specifying the interaction term alone with : in the way you did gives this model matrix with 7 columns:

 model.matrix(~A:B)
  (Intercept) Aa:Ba Ab:Ba Aa:Bb Ab:Bb Aa:Bc Ab:Bc
1           1     1     0     0     0     0     0
2           1     0     1     0     0     0     0
3           1     0     0     1     0     0     0
4           1     0     0     0     1     0     0
5           1     0     0     0     0     1     0
6           1     0     0     0     0     0     1
## attributes omitted from this display

The (Intercept) column for each row is the sum of the subsequent 6 columns, so this model matrix has linearly dependent columns (rank deficiency).

Specifying the standard interaction with *, which gives individual predictor coefficients, results in the correct number of linearly independent columns, 6 (as does specifying all the terms with ~A+B+A:B, not shown).

model.matrix(~A*B)
  (Intercept) Ab Bb Bc Ab:Bb Ab:Bc
1           1  0  0  0     0     0
2           1  1  0  0     0     0
3           1  0  1  0     0     0
4           1  1  1  0     1     0
5           1  0  0  1     0     0
6           1  1  0  1     0     1
## attributes omitted as above

So your attempt to preserve degrees of freedom by omitting "main effects" is what caused a lot of your problem.

People sometimes play with omitting "main effects" or intercepts to make coefficients reported in the initial model summary "easier to interpret." That can pose problems and actually complicate interpretation, as it did for you. For example, some try to remove the intercept with a categorical predictor, to get coefficients that represent mean values rather than an intercept for the reference level and difference for the others. But that only works with a single categorical predictor, as this thread demonstrates.

Separate out the fitting of the model from the interpretation and display of the results, particularly when there are categorical predictors and interactions. All full-rank codings of categorical predictors provide the same fundamental model from which, with additional tools, you can display the results for any desired combination of predictor values.