1
$\begingroup$

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?

$\endgroup$
1
  • $\begingroup$ The way that you specified the interactions with ":" to try to omit "main effects" is one source of your difficulty, as illustrated in an expansion of my answer. $\endgroup$
    – EdM
    Jan 17 at 19:37

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
9
  • $\begingroup$ I see - so it sounds like every random effect in a Mixed model must also be specified as a fixed Main effect. Do you have any insight on how Mean Centering the data might compare with the Deviation Contrasts method in this case? Thank you! $\endgroup$ Jan 16 at 19:05
  • $\begingroup$ @LothartheQuick I'm not sure that "every random effect in a Mixed model must also be specified as a fixed Main effect," but that's how I always do it and I don't know how to interpret coefficients if one doesn't. Mean centering a continuous predictor can make it easier to interpret coefficients of predictors involved in interactions with it, but how do you do that for a categorical predictor like Region? Deviation contrasts are the closest that you can come. I'd suggest obsessing less over the predictor coding and instead learning tools like those in emmeans to illustrate your results. $\endgroup$
    – EdM
    Jan 16 at 19:22
  • $\begingroup$ @LothartheQuick if you are talking about mean-centering the outcome variable then, in a completely balanced design with an ordinary linear regression, you are effectively removing the intercept and thus individual coefficients will represent deviations from the overall mean. But the interaction coefficients will still represent differences from lower-level coefficients, and that approach won't help at all if you need to use a generalized model (as I suspect you will with MRI data). $\endgroup$
    – EdM
    Jan 16 at 19:27
  • $\begingroup$ EdM, Thank you again for the insightful discussion. Yes, I agree that the interpretation of the random effects does not make sense in the absence of a fixed main effect coefficient, but I wasn't sure how the interpretation works within the context of interaction coefficients as of yet. At this point, I'm experimenting with fitting the crossed effects MEM with a limited sample size, so I'm trying to preserve dof's by leaving the main effects terms out if possible. $\endgroup$ Jan 16 at 19:52
  • $\begingroup$ EdM, I apologize for being unclear, and you indeed picked up on the fact that I was implying that I was mean-centering the outcome variables for the regions. I was originally intending to include the Mean "Brain" region as the reference category that would be absorbed into the intercept but for unknown reasons this doesn't work. We only have 3x time points so we've been using a simple linear rather than a generalized model to avoid overfitting the data in this case, and I expect the elimination of the "Brain" region category will help preserve a few more dof's. $\endgroup$ Jan 16 at 20:04

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.