[R-sig-ME] Multivariate mixed effects model

[Afshartous, David]

Thanks for the reference Mareike.  This approach works fine for simple models, but I don't see how to extend to more complex models. A variety of more complex structures can apparently be fit in SAS (Gao et al, Fitting Multivariate Longitudinal Data Using SAS, Paper 187-31).
In the context of the pseudo example below, does anyone  know how to do the following?
## 1) residual error correlated across variables at the same time point?
## 2) residual error correlated across variables at different time point? (e.g., AR1)
## 3) different AR1 coefficient for different variables?

Thanks,
David


## basic bivariate model with correlated intercept random-effects
mv1.lme = lme(y ~ -1 + y1.ind + y2.ind + y1.time + y2.time, random = ~ -1 + y1.ind + y2.ind  | subject,
            data = data.mv, control=nlmeControl(msMaxIter = 500))
## different residual error variance per variable
mv3.lme = lme(y ~ -1 + y1.ind + y2.ind + y1.time + y2.time, random = ~ -1 + y1.ind + y2.ind  | subject,
            data = data.mv, control=nlmeControl(msMaxIter = 500), weights = varIdent(form = ~1 | y1.ind))

## for comparison, basic model in lmer below, but don't think residual variance extensions above currently available:
library("lme4")
mv1.lmer  = lmer(y ~ -1 + y1.ind + y2.ind + y1.time + y2.time +  (0 + y1.ind + y2.ind + y1.time + y2.time | subject) ,
            data = data.mv)

###############
## simulated data:
set.seed(100); n1 = 30; n2 = 4
library("MASS"); library("nlme")
y1 = y2 = matrix(0, n2, n1)
for (i in 1:n1) {
    b.i = mvrnorm(n = 1, mu = c(0,0), Sigma = matrix(c(3, 1, 1, 6), nrow = 2, ncol = 2)) # correlated random effects
   for (j in 1:n2) {
        eps =  mvrnorm(n = 1, mu = c(0,0), Sigma = matrix(c(2, 0, 0, 5), nrow = 2, ncol = 2))  ##  uncorrelated error term
       y1[j,i] = 4 + 2 * j + b.i[1] + eps[1]  ## fixed intercept and sloped effects + random-intercept
       y2[j,i] = 25 + 13 * j + b.i[2] + eps[2] }}
data.mv = data.frame(subject = rep(seq(1,n1), each = n2, 2), time = rep(c(1:n2), n1*2), y = c(as.vector(y1), as.vector(y2)), y1.ind = rep( c(1,0),
    each = n1*n2), y2.ind = rep( c(0,1), each = n1*n2), y1.time = c(rep(seq(1:4), n1), rep(0, n1*n2)), y2.time = c(rep(0, n1*n2), rep(seq(1,4), n1)))



On 6/30/09 4:21 PM, "Mareike Kohlmann" <mareike.kohlmann at stat.uni-muenchen.de> wrote:

Hi,

lme or lmer can be used when the multivariate variables are structured in
the same way as if only a univariate variable was analyzed. A set of dummy
codes need to be created to "flag" the outcomes accordingly.

For more details, see e.g:
Doran, H., Lockwood, J., 2006. Fitting value-added models in R. Journal of
Educational and Behavioral Statistics 31 (2), p. 205-230.

Cheers,
Mareike

***************************************
Mareike Kohlmann
Department of Statistics
Ludwig-Maximilians-University Munich
Germany
E-Mail: mareike.kohlmann at stat.uni-muenchen.de