Can plot interaction means for nlme fit, but not for lme4

Question

Using this dummy dataset:

dummy.data <- data.frame(vaccinated = factor(rep(c("yes", "no"), each = 64)),
  infected = factor(rep(c("yes", "no"), each = 32)),
  animal = factor(rep(1:16, each = 8)),
  tissue = factor(c("blood", "liver", "kidney", "brain")),
  value = runif(128)
  )

This works:

library("nlme")
nlme.model <- as.formula(value ~ vaccinated * infected * tissue)
nlme.fit <- lme(fixed = nlme.model, random = ~1|animal, data = dummy.data)
library("phia")
int.nlme <- interactionMeans(nlme.fit)
plot(int.nlme)

But this doesn't:

library("lme4")  
lmer.model <- as.formula(value ~ vaccinated * infected * tissue + (1 | animal))
lmer.fit <- lmer(formula = lmer.model, data = dummy.data)
library("phia")
int.lmer <- interactionMeans(lmer.fit)
plot(int.lmer)

With the latter, I only get

Error in t.default(M) : argument is not a matrix

from the plot command.

When I look at int.nlme and int.lmer with str, they do look different, but I can't figure out what the problem is. Any input is much appreciated.

Answer 1

The error appears to be generated from phia:::poolse and can be replicated thus:

> phia:::poolse(int.nlme, "adjusted mean","vaccinated")
vaccinated
        no        yes 
0.04483019 0.04483019 
> phia:::poolse(int.lmer, "adjusted mean","vaccinated")
Error in t.default(M) : argument is not a matrix

Still digging...

And I conclude its a bug in the phia package, which has neglected this:

When writing an R package that uses the Matrix package, why do I have to specify Matrix::t() instead of just t()?

as a workaround, and to increase the awesomeness, turn the "covmat" attribute of int.lmer from Matrix classes to standard R matrix classes:

> int.lmer <- interactionMeans(lmer.fit)
> plot(int.lmer)
Error in t.default(M) : argument is not a matrix
> attr(int.lmer, "covmat") = lapply(attr(int.lmer,"covmat"),as.matrix)
> plot(int.lmer)

The plot then works.

Answer 2

I confirm that this is effectively the reason of the bug: the covariance matrix of lmer has to be transformed into a "normal" matrix for the plot to work. Fixing that is priority in my to-do list, and I want to do a series of fixes including that during this summer. But if anyone wants to contribute to accelerate the update, you can suggest a pull request to the repository in Github:

https://github.com/heliosdrm/phia