lme4/lme4 — pull request #144 by runehaubo (2013-10-19 09:18:12+00:00)
Title: Ensuring X is of full rank

I propose adding these changes to automatically remove redundant columns from X to ensure that it has
full rank before fitting the model. It basically comes down to making a LINPACK `qr` decomposition of X and remove columns detected to be linear combinations of other columns down to a tolerance. I have left the orignal check for X rank deficiency, and I have added the control option `ensureXrank=TRUE`, so the automatic check can be disabled if you want to save the time it takes to do a QR on X.

I increased version from 1.1-0 to 1.1-1 to able to distinguish the packages locally, but you might not want that - don't know what your policy is here.

I have updated relevant tests and added some tests related to `anova` on single model fits which are affected.

The main benefit of this change is the ability to fit models such as
````
set.seed(101) ## borrowed from ./inst/tests/test-rank.R
n <- 20
x <- y <- rnorm(n)
z <- rnorm(n)
r <- sample(1:5, size=n, replace=TRUE)
d <- data.frame(x,y,z,r)
d2 <- expand.grid(a=factor(1:4),b=factor(1:4),rep=1:10)
n <- nrow(d2)
d2 <- transform(d2,r=sample(1:5, size=n, replace=TRUE),
z=rnorm(n))
d2 <- subset(d2,!(a=="4" & b=="4"))
fm <- lmer( z ~ a*b + (1|r), data=d2)
design is column rank deficient so dropping 1 coef
anova(fm)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
a 3 3.6087 1.20290 1.3046
b 3 4.4842 1.49473 1.6211
a:b 8 4.6847 0.58559 0.6351
````
As you see, [g]lmer signals that it drops 1 coefficient (and maybe you want an option to turn of that message?). The best thing would naturally be to simulate the `lm` behaviour and have the offending coefficient show up as `NA` in `print` and `summary`, but that involves much deeper changes to the code base:

````
d> fm
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ a * b + (1 | r)
Data: d2
REML criterion at convergence: 406.6967
Random effects:
Groups Name Std.Dev.
r (Intercept) 0.0000
Residual 0.9602
Number of obs: 150, groups: r, 5
Fixed Effects:
(Intercept) a2 a3 a4 b2 b3
-0.29790 -0.15957 0.29289 -0.08277 0.41987 0.61012
b4 a2:b2 a3:b2 a4:b2 a2:b3 a3:b3
0.20171 -0.25529 -0.15419 -0.29456 -0.48281 -0.09069
a4:b3 a2:b4 a3:b4
-0.01160 0.54335 -0.17982
d> lm(z ~ a*b, data=d2)

Call:
lm(formula = z ~ a * b, data = d2)

Coefficients:
(Intercept) a2 a3 a4 b2 b3
-0.29790 -0.15957 0.29289 -0.08277 0.41987 0.61012
b4 a2:b2 a3:b2 a4:b2 a2:b3 a3:b3
0.20171 -0.25529 -0.15419 -0.29456 -0.48281 -0.09069
a4:b3 a2:b4 a3:b4 a4:b4
-0.01160 0.54335 -0.17982 NA
````

Here is another 'real world' example which benefits from the change:

````
d> data(soup, package="ordinal")
d> sureness <- as.numeric(as.character(soup$SURENESS))
d> fm2 <- lmer(sureness ~ PRODID * DAY + (1|RESP), data=soup)
design is column rank deficient so dropping 1 coef
d> anova(fm2)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
PRODID 5 638.06 127.611 41.3472
DAY 1 19.29 19.289 6.2497
PRODID:DAY 4 27.25 6.811 2.2070
````
so now I can get the ANOVA F-test of the interaction.

Hope that you are positive of my proposal.
All the best,
Rune