10

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?

smillig
  • 2,336
  • 28
  • 31
Chris
  • 799
  • 1
  • 7
  • 15

1 Answers1

3

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:

m1 <- lmer(value ~ -1 + variable + variable:gear + variable:carb + (0 + variable | factor(carb))

In that case, you would get the estimates of fixed effects on each variable (for example, variabledrat:gear is the effect of predictor gear on response drat), but you would also get the intercepts for each variable (for example variabledrat for the intercept of response drat) and the random variance of each variable and the correlations between variables:

Groups       Name         Std.Dev. Corr     
 factor(carb) variabledrat 23.80             
              variablempg  24.27    0.20     
              variablehp   23.80    0.00 0.00
 Residual                  23.80

A more detailed description of these methods has been written down by Ben Bolker, as well as the use of MCMCglmmin 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).

Axeman
  • 199
  • 1
  • 11
Charlotte R
  • 438
  • 4
  • 18