non-conformable arguments error from lmer when trying to extract information from the model matrix

I have some longitudinal data from which I'd like to get the predicted means at specified times. The model includes 2 terms, their interaction and a spline term for the time variable. When I try to obtain the predicted means, I get "Error in mm %*% fixef(m4) : non-conformable arguments"

I've used the sleep data set from lmer to illustrate my problem. First, I import the data and create a variable "age" for my interaction

sleep <- as.data.frame(sleepstudy)  #get the sleep data
# create fake variable for age with 3 levels
set.seed(1234567)
sleep$age <- as.factor(sample(1:3,length(sleep),rep=TRUE))

Then I run my lmer model

library(lme4)
library(splines)
m4 <- lmer(Reaction ~ Days + ns(Days, df=4) + age + Days:age + (Days | Subject), sleep)

Finally, I create the data and matrix needed to obtain predicted means

#new data frame for predicted means
d <- c(0:9)  # make a vector of days = 0 to 9 to obtain predictions for each day
newdat <- as.data.frame(cbind(Days=d, age=rep(c(1:3),length(d))))
newdat$Days <- as.numeric(as.character(newdat$Days))
newdat$age <- as.factor(newdat$age)

# create a matrix 
mm<-model.matrix(~Days + ns(Days, df=4) + age + Days:age, newdat)  
newdat$pred<-mm%*%fixef(m4)

It's at this point that I get the error: Error in mm %*% fixef(m4) : non-conformable arguments

I can use predict to get the means

newdat$pred <- predict(m4, newdata=newdat, re.form=NA)

which works fine, but I want to be able to calculate a confidence interval, so I need a conformable matrix.

I read somewhere that the problem may be that lmer creates aliases (I can't find that post). This comment was made with regards to not being able to use effect() for a similar task. I couldn't quite understand how to overcome this problem. Moreover, I recall that post was a little old and hoped the alias problem may no longer be relevant.

If anyone has a suggestion for what I may be doing wrong, I'd appreciate the feedback. Thanks.

Solution

There are a couple of things here.

you need to drop columns to make your model matrix commensurate with the fixed effect vector that was actually fitted (i.e., commensurate with the model matrix that was actually used for fitting, after dropping collinear columns)
for additional confusion, you happened to only sample ages 2 and 3 (out of a possible {1,2,3})

I've cleaned up the code a little bit ...

library("lme4")
library("splines")
sleep <- sleepstudy  #get the sleep data
set.seed(1234567)
## next line happens to sample only 2 and 3 ...
sleep$age <- as.factor(sample(1:3,length(sleep),rep=TRUE))
length(levels(sleep$age))  ## 2

Fit model:

m4 <- lmer(Reaction ~ Days + ns(Days, df=4) +
    age + Days:age + (Days | Subject), sleep)
## message; fixed-effect model matrix is 
##    rank deficient so dropping 1 column / coefficient

Check fixed effects:

f1 <- fixef(m4)
length(f1)  ## 7
f2 <- fixef(m4,add.dropped=TRUE)
length(f2)  ## 8

We could use this extended version of the fixed effects (which has an NA value in it), but this would just mess us up by propagating NA values through the computation ...

Check model matrix:

X <- getME(m4,"X")
ncol(X)  ## 7
(which.dropped <- attr(getME(m4,"X"),"col.dropped"))
## ns(Days, df = 4)4 
##             6

New data frame for predicted means

d <- 0:9  
## best to use data.frame() directly, avoid cbind()
##   generate age based on *actual* levels in data
newdat <- data.frame(Days=d,
   age=factor(rep(levels(sleep$age),length(d))))

Create a matrix:

mm <- model.matrix(formula(m4,fixed.only=TRUE)[-2], newdat)
mm <- mm[,-which.dropped]   ## drop redundant columns
## newdat$pred <- mm%*%fixef(m4)    ## works now

Added by sianagh: Code to obtain confidence intervals and plot the data:

predFun <- function(x) predict(x,newdata=newdat,re.form=NA)
newdat$pred <- predFun(m4)
bb <- bootMer(m4,
   FUN=predFun,
    nsim=200)  
## nb. this produces an error message on its first run, 
## but not on subsequent runs (using the development version of lme4)
bb_ci <- as.data.frame(t(apply(bb$t,2,quantile,c(0.025,0.975))))
names(bb_ci) <- c("lwr","upr")
newdat <- cbind(newdat,bb_ci)

Plot:

plot(Reaction~Days,sleep)
with(newdat,
    matlines(Days,cbind(pred,lwr,upr),
            col=c("red","green","green"),
            lty=2,
            lwd=c(3,2,2)))