Hard to tell for sure without a reproducible example: How to make a great R reproducible example?
But, guessing: these sorts of problems are generally due to collinearity in the design matrix. Centering your continuous predictor (intdiff
) may help. You can also explore the design matrix directly
X <- model.matrix( ~ stress_limit * word_position * follows, data)
Collinearity between pairs: cor(X)
. Unfortunately I don't have a suggestion for detecting multi-collinearity (i.e. not between pairs, but between combinations of >2 predictors) off the top of my head, although you can look into the tools for computing variance inflation factors (e.g. library("sos"); findFn("VIF")
).
As a cross-check, lme
should also be able to handle your model:
library(nlme)
lme(intdiff ~ stress_limit * word_position * follows,
random=~1|speaker, data=data)
When I run your test data in the development version of lme4 (available on github) I get Error in lmer(intdiff ~ stress_limit * word_position * follows + (1 | : rank of X = 5 < ncol(X) = 12
. On the other hand, with this small an input data set (6 observations), there's no possible way you could fit 12 parameters. It's a little harder to tell exactly where your problem is. Do all 12 combinations of your 3 variables actually occur in your data? If some are missing, then you need to follow the advice given in the development version's help:
Unlike some simpler modeling frameworks such as ‘lm’ and ‘glm’ which automatically detect perfectly collinear predictor variables, ‘[gn]lmer’ cannot handle design matrices of less than full rank. For example, in cases of models with interactions that have unobserved combinations of levels, it is up to the user to define a new variable (for example creating ‘ab’ within the data from the results of ‘droplevels(interaction(a,b))’).
In particular, you can fit this model as follows:
data <- transform(data,
allcomb=interaction(stress_limit,word_position,follow,drop=TRUE))
lme(intdiff ~ allcomb, random=~1|speaker, data=data)
This will give you a one-way ANOVA treating the unique combinations of levels that are actually present in the data as the categories. You'll have to figure out for yourself what they mean.
The alternative is to reduce the number of interactions in the model until you get to a set that don't have any missing combinations; if you're lucky (stress_limit+word_position+follow)^2
(all two-way interactions) will work, but you might have to reduce the model still farther (e.g. stress_limit + word_position*follow
).
Another way to test this is to use lm()
on your proposed models and check that there are no NA
values in the estimated coefficients.
The main thing you will be losing in these ways is convenience/ease of interpretation, because the parameters for the missing combinations couldn't have been estimated from the data anyway ...