Correlations Fixed Effect

Question

I am fitting a generalised mixed linear model using the lme4 package and glmer. Each time I fit the model, I get this error message before the model even starts fitting:

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

Surprisingly, the question of WHICH effect gets dropped depends on the order of the variables in the formula. If I use y ~ a + b + c + a*b + a*c then, a*c gets dropped, if I do y ~ a + b + c + a*c + a*b, then a*b gets dropped. However, cor(a*c, a*b)=0.03, which contradicts my intuition of multicollinearity. So I started investigating what's going on.

During this investigation I noted the following. If I use summary(model) to find the correlations of the fixed effects, the correlations are MUCH higher than when I just compute the correlation matrix using cor(x).

Q1: Why are these two tables different?

Q2: summary(model) only shows the variables which are NOT dropped. Can I recreate the table including effects that have been dropped?

Answer 1

Your assumption of multi-collinearity is correct, but it's not necessary for multi-collinearity that a*c and a*b are correlated in themselves.

Multi-collinearity occurs if any linear combination of terms can completely determine any of the others.

For example, suppose I have a dataset with independent predictors $X_1\sim N(0,1)$ and $X_2\sim N(0,1)$. Then calculate $X_3=X_1+X_2$ and $X_4=X_1-X_2$.

Now $X_3$ and $X_4$ are uncorrelated, but given $X_1$ and $X_3$ I can completely determine $X_4$ (because $X_4=2X_1-X_3$).

So I cannot estimate a regression with $X_1$, $X_3$ and $X_4$ as predictors. Intuitively I hope it's obvious that you wouldn't learn anything new about an outcome by supplying $X_4$ given $X_1$ and $X_3$, likewise if you already have $X_1$ and $X_4$ you don't learn anything by adding $X_3$.

Also from a mathematical point of view you can't invert $X^TX$ (where $X$ is the design matrix whose columns include $X_1$, $X_3$, and $X_4$).

So to answer your second question, it doesn't make sense to try to estimate the coefficients for the dropped effects, when the other effects are also in the model.

To see it using R:

# choose a sample size
N=1000

# make independent predictors
> df <- data.frame(x1=rnorm(N), x2=rnorm(N))

# make linearly dependent predictors
> df$x3 <- df$x1 - df$x2
> df$x4 <- df$x1 + df$x2

# make an outcome
> df$y <- rnorm(N)

# check that x3 and x4 are not correlated
> cor(df$x3, df$x4)
[1] 0.05538154

# try to estimate the effect of x3 and x4 on y, given x1.
> lm(data=df , y~x1+x3+x4)

Call:
lm(formula = y ~ x1 + x3 + x4, data = df)

Coefficients:
(Intercept)           x1           x3           x4  
  -0.009205     0.007510    -0.016671           NA  

> lm(data=df , y~x1+x4+x3)

Call:
lm(formula = y ~ x1 + x4 + x3, data = df)

Coefficients:
(Intercept)           x1           x4           x3  
  -0.009205    -0.025833     0.016671           NA  

> lm(data=df , y~x4+x3+x1)

Call:
lm(formula = y ~ x4 + x3 + x1, data = df)

Coefficients:
(Intercept)           x4           x3           x1  
  -0.009205     0.003755    -0.012916           NA  

# Check the model matrix:
> M=model.matrix(data=df , y~x1+x4+x3)
> solve(t(M)%*%M)
Error in solve.default(t(M) %*% M) : 
  system is computationally singular: reciprocal condition number = 1.09216e-16

R needs to drop a term each time before it can run the regression. It doesn't matter which term it drops, so it picks the last one. The model is statistically identical in each case.