I'm trying to fit a multivariate (i.e., multiple response) mixed model in R
. Aside from the ASReml-r
and SabreR
packages (which require external software), it seems this is only possible in MCMCglmm
. In the paper that accompanies the MCMCglmm
package (pp.6), Jarrod Hadfield describes the process of fitting such a model as like reshaping multiple response variables into one long-format variable and then suppressing the overall intercept. My understanding is that suppressing the intercept changes the interpretation of the coefficient for each level of the response variable to be the mean for that level. Given the above, is it therefore possible to fit a multivariate mixed model using lme4
? For example:
data(mtcars)
library(reshape2)
mtcars <- melt(mtcars, measure.vars = c("drat", "mpg", "hp"))
library(lme4)
m1 <- lmer(value ~ -1 + variable:gear + variable:carb + (1 | factor(carb)),
data = mtcars)
summary(m1)
# Linear mixed model fit by REML
# Formula: value ~ -1 + variable:gear + variable:carb + (1 | factor(carb))
# Data: mtcars
# AIC BIC logLik deviance REMLdev
# 913 933.5 -448.5 920.2 897
# Random effects:
# Groups Name Variance Std.Dev.
# factor(carb) (Intercept) 509.89 22.581
# Residual 796.21 28.217
# Number of obs: 96, groups: factor(carb), 6
#
# Fixed effects:
# Estimate Std. Error t value
# variabledrat:gear -7.6411 4.4054 -1.734
# variablempg:gear -1.2401 4.4054 -0.281
# variablehp:gear 0.7485 4.4054 0.170
# variabledrat:carb 5.9783 4.7333 1.263
# variablempg:carb 3.3779 4.7333 0.714
# variablehp:carb 43.6594 4.7333 9.224
How would one interpret the coefficients in this model? Would this method also work for generalized linear mixed models?
Best Answer
Your idea is good, but in your example, you forgot to model different intercept and different random variances for each trait, so your output is not interpretable as is. A correct model would be:
In that case, you would get the estimates of fixed effects on each variable (for example,
variabledrat:gear
is the effect of predictorgear
on responsedrat
), but you would also get the intercepts for each variable (for examplevariabledrat
for the intercept of responsedrat
) and the random variance of each variable and the correlations between variables:A more detailed description of these methods has been written down by Ben Bolker, as well as the use of
MCMCglmm
in a Bayesian framework. Another new package,mcglm
can also handle multivariate models, even with non-normal responses, but you have to code your random design matrices. A tutorial should be available soon (see the R help page).