15

I wonder how to fit multivariate linear mixed model with lme4. I fitted univariate linear mixed models with the following code:

library(lme4)
lmer.m1 <- lmer(Y1~A*B+(1|Block)+(1|Block:A), data=Data)
summary(lmer.m1)
anova(lmer.m1)

lmer.m2 <- lmer(Y2~A*B+(1|Block)+(1|Block:A), data=Data)
summary(lmer.m2)
anova(lmer.m2)

I'd like to know how to fit multivariate linear mixed model with lme4. The data is below:

Block A B    Y1    Y2
 1    1 1 135.8 121.6
 1    1 2 149.4 142.5
 1    1 3 155.4 145.0
 1    2 1 105.9 106.6
 1    2 2 112.9 119.2
 1    2 3 121.6 126.7
 2    1 1 121.9 133.5
 2    1 2 136.5 146.1
 2    1 3 145.8 154.0
 2    2 1 102.1 116.0
 2    2 2 112.0 121.3
 2    2 3 114.6 137.3
 3    1 1 133.4 132.4
 3    1 2 139.1 141.8
 3    1 3 157.3 156.1
 3    2 1 101.2  89.0
 3    2 2 109.8 104.6
 3    2 3 111.0 107.7
 4    1 1 124.9 133.4
 4    1 2 140.3 147.7
 4    1 3 147.1 157.7
 4    2 1 110.5  99.1
 4    2 2 117.7 100.9
 4    2 3 129.5 116.2

Thank in advance for your time and cooperation.

MYaseen208
  • 22,666
  • 37
  • 165
  • 309
  • 2
    It may be possible to do this by 'melting' the data set (i.e. making Y1 and Y2 separate observations with a common 'ID' variable) and then fitting a model with ID as a random effect. Don't have time to elaborate now. You might want to ask this on the r-sig-mixed-models list. – Ben Bolker Oct 20 '11 at 13:24
  • @Ben Bolker: Thanks for your comment. I'm waiting for your reply. Thanks – MYaseen208 Oct 25 '11 at 03:18
  • 2
    I've put up some more details at http://rpubs.com/bbolker/3336 – Ben Bolker Jan 05 '13 at 20:17

3 Answers3

14

This can sometimes be faked satisfactorily in nlme/lme4 by simply reformatting your data like

require(reshape)
Data = melt(data, id.vars=1:3, variable_name='Y')
Data$Y = factor(gsub('Y(.+)', '\\1', Data$Y))
> Data
  Block A B Y value
1     1 1 1 1 135.8
2     1 1 2 1 149.4
3     1 1 3 1 155.4
4     1 2 1 1 105.9
5     1 2 2 1 112.9
6     1 2 3 1 121.6
...

and then including the new variable Y in your linear mixed model.

However, for true Multivariate Generalized Linear Mixed Models (MGLMM), you will probably need the sabreR package or similar. There is also an entire book to accompany the package, Multivariate Generalized Linear Mixed Models Using R. If you have a proxy to a subscribing institution, you might even be able to download it for free from http://www.crcnetbase.com/isbn/9781439813270. I would refer you there for any further advice, as this is a meaty topic and I am very much a novice.

John Colby
  • 22,169
  • 4
  • 57
  • 69
  • Thanks for your reply. I'd highly appreciate if you give more guidelines to fit this kind of model with lme4 or nlme. Thanks – MYaseen208 Oct 28 '11 at 06:03
  • You have the basic syntax, but determining what specific model to fit will depend completely on your data and the structure of your experiment. This probably isn't the right forum for that. If you aren't familiar with mixed models, *Mixed Effects Models in S and S-Plus* is the book that accompanies the `nlme` package, and might be a good place to start. – John Colby Oct 28 '11 at 15:14
5

lmer and its elder sibling lme are inherently "one parameter left of ~". Have a look at the car packages; it offers no off-the shelf repeated measurement support, but you will find a few comments on the subject by searching the R list:

John Fox on car package

IRTFM
  • 258,963
  • 21
  • 364
  • 487
Dieter Menne
  • 10,076
  • 44
  • 67
1

@John's answer above should be largely right. You add a dummy variable (ie--the factor variable Y) to the model. Here you have 3 subscripts i= 1...N for observations, j=1,...,4 for blocks, and h=1,2 for the dependent var. But you also need to force the level 1 error term to 0 (or to near zero), which I'm not sure lme4 does. Ben Bolker might provide more information. This is described more in Goldstein (2011) Chap 6 and Chap 7 for latent multivariate models.

IE

Y_hij = \beta_{01} z_{1ij} + \beta_{02} z_{2ij} + \beta X + u_{1j} z_{1ij} + u_{2j} z_{2ij}

So:

require(reshape2)
Data = melt(data, id.vars=1:3, variable_name='Y')
Data$Y = factor(gsub('Y(.+)', '\\1', Data$Y))

m1 <- lmer(value ~ Y + A*B + (1|Block) + (1|Block*A), data= Data)
# not sure how to set the level 1 variance to 0, @BenBolker
# also unclear to me if you're requesting Y*A*B instead of Y + A*B
alexwhitworth
  • 4,839
  • 5
  • 32
  • 59