0
$\begingroup$

If I fit an lmer() model with correlated fixed effects, by default it gives me a fit warning fixed-effect model matrix is rank deficient and returns the model with one or more parameters dropped.

For example:

set.seed(1030)
d = data.frame(y=rnorm(100), x1=c('a', 'b'), x2=c('c', 'd'), k=1:5)
lmer(y ~ x1 + x2 + (1 + x1 + x2 | k), d)

(plus convergence warnings etc)

This fits a model with one fixed-effect parameter dropped. Is it possible to also force it to drop the slopes that correspond to the dropped column(s)?

I found instructions on how to identify the bad columns manually, but my situation is this: I have an interaction of multilevel factors, where some combinations of factor levels never co-occur in the data. The factors are sum-coded, and I'd like to just fit the maximal model without hand-coding new columns for each contrast parameter, because I'm using lsmeans to calculate predicted marginal means. If I can keep the original factor coding then I can get also get an easily-interpretable output from lsmeans.

Edit: an example like my data

This is what I'm doing now:

set.seed(1030)

# level 3 of 'a' only occurs with level 3 of 'b'
x = expand.grid(a=1:3, b=1:3)
x = subset(x, !(a == 3 & b != 3))

d = data.frame(y=rnorm(700),
               a=factor(x$a),
               b=factor(x$b),
               k=factor(1:5))

library(lme4)
library(lsmeans)

# tmp for speed
lsm.options(disable.pbkrtest=TRUE)

# this model includes random slopes for factor level combinations that don't exist
m1 = lmer(y ~ a * b + (1 + a * b | k), data=d)

lsmeans(m1, ~ a | b)

# and then pairs() for comparisons

Is this (below) what you're suggesting, or is this not quite right?

# this model has a different random parameterization, but 
# doesn't include slopes for non-existent combinations
m2 = lmer(y ~ a * b + (1 + droplevels(interaction(a, b)) | k), data=d)

lsmeans(m2, ~ a | b)
# produces the same results as above
$\endgroup$
  • $\begingroup$ The results from lsmeans won't depend on how factors are coded. But some of the results may not be estimable -- NAs will be returned for those, again regardless of how they are coded. $\endgroup$ – rvl Oct 30 '15 at 21:59
  • $\begingroup$ Right, lsmeans is giving me NAs exactly where I expect it to for the non-estimable cells & contrasts. I just mean that I want to be able to use lsmeans(model, ~ factor), where factor is a 3-level sum-coded factor that both lmer and lsmeans know how to deal with, rather than lsmeans(model, ~ factor_1 + factor_2), where factor_1 is a column coded as 1 for level 1, 0 for level 2, and -1 for level 3; and factor_2 is a column coded as 0 for level 1, 1 for level 2, and -1 for level 3. $\endgroup$ – user2034412 Oct 31 '15 at 1:00
  • $\begingroup$ If I need to make those columns by hand so that I can exclude some combinations of them for the interaction slopes, then lsmeans won't know which combinations of values correspond to which factor levels, and I'll have to interpret that myself. $\endgroup$ – user2034412 Oct 31 '15 at 1:02
  • $\begingroup$ Why not just create a new factor using interaction(factor1, factor2)? Or am I missing something? $\endgroup$ – rvl Oct 31 '15 at 2:03
  • $\begingroup$ I'm not totally sure what you mean - I edited my question to include an example like my data. Thank you for your help, by the way! I really appreciate it. $\endgroup$ – user2034412 Oct 31 '15 at 18:17
2
$\begingroup$

I suggest something like this:

newd = transform(d, ab = interaction(a, b))
m2 = lmer(y ~ ab + (1 + ab | k), data = newd)
lsm = lsmeans(m2, ~ ab)
contrast(lsm, list(c1 = c(...), c2 = c(...), ...))

In other words, use the new factor ab, which consists of the combinations of a and b that actually occur in your dataset, as the one factor to use in the model in place of every usage of a and b. In the contrast call, put the coefficients for meaningful comparisons or contrasts among the levels of ab.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.