What is rank deficiency, and how to deal with it?

Question

Fitting a logistic regression using lme4 ends with

Error in mer_finalize(ans) : Downdated X'X is not positive definite. 

A likely cause of this error is apparently rank deficiency. What is rank deficiency, and how should I address it?

Answer 1

Rank deficiency in this context says there is insufficient information contained in your data to estimate the model you desire. It stems from many origins. I'll talk here about modeling in a fairly general context, rather than explicitly logistic regression, but everything still applies to the specific context.

The deficiency may stem from simply too little data. In general, you cannot uniquely estimate n parameters with less than n data points. That does not mean that all you need are n points, as if there is any noise in the process, you would get rather poor results. You need more data to help the algorithm to choose a solution that will represent all of the data, in a minimum error sense. This is why we use least squares tools. How much data do you need? I was always asked that question in a past life, and the answer was more than you have, or as much as you can get. :)

Sometimes you may have more data than you need, but some (too many) points are replicates. Replication is GOOD in the sense that it helps to reduce the noise, but it does not help to increase numerical rank. Thus, suppose you have only two data points. You cannot estimate a unique quadratic model through the points. A million replicates of each point will still not allow you to fit more than a straight line, through what are still only effectively a pair of points. Essentially, replication does not add information content. All it does is decrease noise at locations where you already have information.

Sometimes you have information in the wrong places. For example, you cannot fit a two dimensional quadratic model if all you have are points that all lie in a straight line in two dimensions. That is, suppose you have points scattered only along the line x = y in the plane, and you wish to fit a model for the surface z(x,y). Even with zillions of points (not even replicates) will you have sufficient information to intelligently estimate more than a constant model. Amazingly, this is a common problem that I've seen in sampled data. The user wonders why they cannot build a good model. The problem is built into the very data they have sampled.

Sometimes it is simply choice of model. This can be viewed as "not enough data", but from the other side. You wish to estimate a complicated model, but have provided insufficient data to do so.

In all of the above instances the answer is to get more data, sampled intelligently from places that will provide information about the process that you currently lack. Design of experiments is a good place to start.

However, even good data is sometimes inadequate, at least numerically so. (Why do bad things happen to good data?) The problem here may be model related. It may lie in nothing more than a poor choice of units. It may stem from the computer programming done to solve the problem. (Ugh! Where to start?)

First, lets talk about units and scaling. Suppose I try to solve a problem where one variable is MANY orders of magnitude larger than another. For example, suppose I have a problem that involves my height and my shoe size. I'll measure my height in nanometers. So my height would be roughly 1.78 billion (1.78e9) nanometers. Of course, I'll choose to measure my shoe size in kilo-parsecs, so 9.14e-21 kilo-parsecs. When you do regression modeling, linear regression is all about linear algebra, which involves linear combinations of variables. The problem here is these numbers are different by hugely many orders of magnitude (and not even the same units.) The mathematics will fail when a computer program tries to add and subtract numbers that vary by so many orders of magnitude (for a double precision number, that absolute limit is roughly 16 powers of 10.)

The trick is usually to use common units, but on some problems even that is an issue when variables vary by too many orders of magnitude. More important is to scale your numbers to be similar in magnitude.

Next, you may see problems with big numbers and small variation in those numbers. Thus, suppose you try to build a moderately high order polynomial model with data where your inputs all lie in the interval [1,2]. Squaring, cubing, etc., numbers that are on the order of 1 or 2 will cause no problems when working in double precision arithmetic. Alternatively, add 1e12 to every number. In theory, the mathematics will allow this. All it does is shift any polynomial model we build on the x-axis. It would have exactly the same shape, but be translated by 1e12 to the right. In practice, the linear algebra will fail miserably due to rank deficiency problems. You have done nothing but translate the data, but suddenly you start to see singular matrices popping up.

Usually the comment made will be a suggestion to "center and scale your data". Effectively this says to shift and scale the data so that it has a mean near zero and a standard deviation that is roughly 1. That will greatly improve the conditioning of most polynomial models, reducing the rank deficiency issues.

Other reasons for rank deficiency exist. In some cases it is built directly into the model. For example, suppose I provide the derivative of a function, can I uniquely infer the function itself? Of course not, as integration involves a constant of integration, an unknown parameter that is generally inferred by knowledge of the value of the function at some point. In fact, this sometimes arises in estimation problems too, where the singularity of a system is derived from the fundamental nature of the system under study.

I surely left out a few of the many reasons for rank deficiency in a linear system, and I've prattled along for too long now. Hopefully I managed to explain those I covered in simple terms, and a way to alleviate the problem.

Answer 2

For the definition of the rank of a matrix, you can refer to any good textbook on linear algebra, or have a look at the Wikipedia page.

A $n \times p$ matrix $X$ is said to be full rank if $n \geq p$, and its columns are not a linear combination of each other. In that case, the $p \times p$ matrix $X^TX$ is positive definite, which implies that it has an inverse $(X^TX)^{-1}$.

If $X$ is not full rank, one of the columns is fully explained by the others, in the sense that it is a linear combination of the others. A trivial example is when a column is duplicated. This can also happen if you have a 0-1 variable and a column consists of only 0 or only 1. In that case, the rank of the matrix $X$ is less than $n$ and $X^TX$ has no inverse.

Since the solution of many regression problems (including logistic regression) involves the computation intermediate $(X^TX)^{-1}$, it is then impossible to estimate the parameters of the model. Out of curiosity, you can check here how this term is involved in the formula of multiple linear regression.

That was it for absolute rank deficiency. But sometimes the problem shows up when the matrix $X$ is "almost" not full rank, as extensively detailed by @woodchips. This problem is usually referred to as multicollinearity. This issue is fairly common, you can find out more on how to deal with it on related posts here and there.

Answer 3

user974's answer is fantastic from a modelling perspective and gui11aume's answer is fantastic from a mathematical perspective. I want to refine the former answer strictly from a mixed modelling perspective: specifically a generalized mixed modelling (GLMM) perspective. As you can see, you have referenced the R function mer_finalize which is in the fantastic lme4 package. You also say that you are fitting a logistic regression model.

There are many issues that crop up with such types of numerical algorithms. The issue of the matrix structure of the model matrix of fixed effects is certainly worth considering, as user974 alluded to. But this is very easy to assess, simply calculated the model.matrix of your formula= and data= arguments in a model and take its determinant using the det function. Random effects, however, greatly complicate the interpretation, numerical estimation routine, and inference on the fixed effects (what you typically think of as regression coefficients in a "regular" regression model).

Suppose in the simplest case you are only fitted a random intercepts model. Then you are basically considering there to be thousands of unmeasured sources of heterogeneity that are held constant in repeated measures within clusters. You estimate a "grand" intercept, but account for the heterogeneity by assuming that the cluster-specific intercepts have some mean0 normal distribution. The intercepts are iteratively estimated and used to update model effects until convergence is achieved (the log likelihood--or an approximation of it--is maximized). The mixed model is very easy to envision, but mathematically the likelihood is very complex and prone to issues with singularities, local minimae, and boundary points (odds ratios = 0 or infinity). Mixed models do not have quadratic likelihoods like canonical GLMs.

Unfortunately, Venerables and Ripley did not invest much into diagnostics for converge-fails like yours. It is practically impossible to even speculate at the myriad of possible errors leading to such a message. Consider, then, the types of diagnostics I use below:

1. How many observations are there per cluster?
2. What are the results from a marginal model fit using GEE?
3. What is the ICC of the clusters? Is the within cluster heterogeneity close to the between cluster heterogeneity?
4. Fit a 1-step estimator and look at the estimated random effects. Are they approximately normal?
5. Fit the Bayesian mixed model and look at the posterior distribution for the fixed effects. Do they appear to have an approximately normal distribution?
6. Look at a panel plot of the clusters showing exposure or regressor of interest against the outcome using a smoother. Are the trends consistent and clear or are there many possible ways that such trends may be explained? (e.g. what is the "risk" among unexposed subjects, does the exposure appear protective or harmful?)
7. Is it possible to restrict the sample to subjects only having a sufficient number of observations-per-person (say, n=5 or n=10) to estimate the "ideal" treatment effect?

Alternately, you can consider some different modelling approaches:

1. Are there few enough clusters or time points that you can use a fixed effect (such as group indicators or a polynomial time effect) to model the cluster level/autoregressive heterogeneity?
2. Is a marginal model appropriate (using a GEE to improve efficiency of standard error estimation, but still using only fixed effects)
3. Could a Bayesian model with an informative prior on the random effects improve estimation?

Answer 4

I also kept getting warnings about rank deficiency. In my case it seemed to be caused by the fact that I had a categorical variable where some of the categories were empty/not represented in the training set. When I created a category called other for the different low-frequency categories, I got rid of the warning.

This was when building a linear regression with the train() function from caret. Here is the warning:

Warning messages:
In predict.lm(modelFit, newdata) :
prediction from a rank-deficient fit may be misleading