Lme4::Lmer Reports "Fixed-Effect Model Matrix Is Rank Deficient", Do I Need a Fix and How To

lme4::lmer reports fixed-effect model matrix is rank deficient, do I need a fix and how to?

You are slightly over-concerned with the warning message:

fixed-effect model matrix is rank deficient so dropping 7 columns / coefficients.

It is a warning not an error. There is neither misuse of lmer nor ill-specification of model formula, thus you will obtain an estimated model. But to answer your question, I shall strive to explain it.

During execution of lmer, your model formula is broken into a fixed effect formula and a random effect formula, and for each a model matrix is constructed. Construction for the fixed one is via the standard model matrix constructor model.matrix; construction for the random one is complicated but not related to your question, so I just skip it.

For your model, you can check what the fixed effect model matrix looks like by:

fix.formula <- F2_difference ~ sex + nasal + type + vowelLabel + 
               type * vowelLabel + nasal * type

X <- model.matrix (fix.formula, data.df)

All your variables are factor so X will be binary. Though model.matrix applies contrasts for each factor and their interaction, it is still possible that X does not end up with full column rank, as a column may be a linear combination of some others (which can either be precise or numerically close). In your case, some levels of one factor may be nested in some levels of another.

Rank deficiency can arise in many different ways. The other answer shares a CrossValidated answer offering substantial discussions, on which I will make some comments.

• For case 1, people can actually do a feature selection model via say, LASSO.
• Cases 2 and 3 are related to the data collection process. A good design of experiment is the best way to prevent rank-deficiency, but for many people who build models, the data are already there and no improvement (like getting more data) is possible. However, I would like to stress that even for a dataset without rank-deficiency, we can still get this problem if we don't use it carefully. For example, cross-validation is a good method for model comparison. To do this we need to split the complete dataset into a training one and a test one, but without care we may get a rank-deficient model from the training dataset.
• Case 4 is a big problem that could be completely out of our control. Perhaps a natural choice is to reduce model complexity, but an alternative is to try penalized regression.
• Case 5 is a numerical concern leading to numerical rank-deficiency and this is a good example.
• Cases 6 and 7 tell the fact that numerical computations are performed in finite precision. Usually these won't be an issue if case 5 is dealt with properly.

So, sometimes we can workaround the deficiency but it is not always possible to achieve this. Thus, any well-written model fitting routine, like lm, glm, mgcv::gam, will apply QR decomposition for X to only use its full-rank subspace, i.e., a maximum subset of X's columns that gives a full-rank space, for estimation, fixing coefficients associated with the rest of the columns at 0 or NA. The warning you got just implies this. There are originally ncol(X) coefficients to estimate, but due to deficiency, only ncol(X) - 7 will be estimated, with the rest being 0 or NA. Such numerical workaround ensures that a least squares solution can be obtained in the most stable manner.

To better digest this issue, you can use lm to fit a linear model with fix.formula.

fix.fit <- lm(fix.formula, data.df, method = "qr", singular.ok = TRUE)

method = "qr" and singular.ok = TRUE are default, so actually we don't need to set it. But if we specify singular.ok = FALSE, lm will stop and complain about rank-deficiency.

lm(fix.formula, data.df, method = "qr", singular.ok = FALSE)
#Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
#  singular fit encountered

You can then check the returned values in fix.fit.

p <- length(coef)
coef <- fix.fit$coef
no.NA <- sum(is.na(coef))
rank <- fix.fit$rank

It is guaranteed that p = ncol(X), but you should see no.NA = 7 and rank + no.NA = p.

Exactly the same thing happens inside lmer. lm will not report deficiency while lmer does. This is in fact informative, as too often, I see people asking why lm returns NA for some coefficients.

Update 1 (2016-05-07):

Let me see if I have this right: The short version is that one of my predictor variables is correlated with another, but I shouldn't worry about it. It is appropriate to use factors, correct? And I can still compare models with anova or by looking at the BIC?

Don't worry about the use of summary or anova. Methods are written so that the correct number of parameters (degree of freedom) will be used to produce valid summary statistics.

Update 2 (2016-11-06):

Let's also hear what package author of lme4 would say: rank deficiency warning mixed model lmer. Ben Bolker has mentioned caret::findLinearCombos, too, particularly because the OP there want to address deficiency issue himself.

Update 3 (2018-07-27):

Rank-deficiency is not a problem for valid model estimation and comparison, but could be a hazard in prediction. I recently composed a detailed answer with simulated examples on CrossValidated: R lm, Could anyone give me an example of the misleading case on “prediction from a rank-deficient”? So, yes, in theory we should avoid rank-deficient estimation. But in reality, there is no so-called "true model": we try to learn it from data. We can never compare an estimated model to "truth"; the best bet is to choose the best one from a number of models we've built. So if the "best" model ends up rank-deficient, we can be skeptical about it but probably there is nothing we can do immediately.

rank deficiency warning mixed model lmer

This is not really specifically a linear-mixed model problem.

It comes down to the fact that you can't estimate an interaction if you don't have any treatment happening in the 'before' period (year 0).

Simplest possible example:

(dd <- data.frame(y=1:3,treat=c(0,0,1),year=c(0,1,1)))

##   y treat year
## 1 1     0    0
## 2 2     0    1
## 3 3     1    1

Fit the model:

lm(y~treat*year,dd) ## == year+treat+year:treat
## Call:
## lm(formula = y ~ treat * year, data = dd)
## 
## Coefficients:
## (Intercept)        treat         year   treat:year  
##           1            1            1           NA  

lm doesn't warn you, but it effectively does the same thing as lmer by removing the extra, collinear column and giving its parameter an NA value. If you try caret::findLinearCombos(dd[c("year","treat")]) you won't get anything back (year and treat are not perfectly collinear), but if you look at the model matrix that R constructs to include the treatment column, you will get something:

X <- model.matrix(~year*treat,dd)
caret::findLinearCombos(X)
## $linearCombos
## $linearCombos[[1]]
## [1] 4 3
## $remove
## [1] 4

This experimental design simply doesn't allow you to estimate the interaction. If you remove it from the formula (use year+treat instead of year*treat) you'll get the same answer, but without the message. Alternatively, in a typical "before-after-control-impact" design (in environmental impact assessment), you would label the individuals who would be getting the treatment as "impact" or "treated" individuals even in year 0; then the interaction would be your actual estimated effect of treatment.

Related Topics