[R-sig-ME] Multivariate mixed effects model

[lorenz.gygax at art.admin.ch]

Dear David and others interested,

I am also trying to get to terms with mixed-effects models having several symultaneous outcomes. Thus, I was very pleased to see the Reference that Mareike provided (and slightly disappointed when I saw how little of the publication actually referred to the multiple outcomes - it nevertheless seems/seemed to be a good starting point).

When running a couple of toy examples the following points struck me. I am not quite sure which of those are actually relevant and which ones just seemed odd to me because my understanding of the corresponding math does not reach far enough.

> This approach works fine for simple models, but I don't see how to
> extend to more complex models.

I agree (see below). And my hunch is that we are asking too much and that lme/lmer can not currently be "tricked" into including all aspects of interest to us. This is just not the application, these methods were programmed to deal with.

> A variety of more complex structures can apparently be fit in SAS
> (Gao et al, Fitting Multivariate Longitudinal Data Using SAS, Paper
> 187-31).

As complex as what you are aiming at? What we would like seems way more complex than the standard mixed models ... which explains why the approach in lme/lmer is somewhat limited (or perhaps I am missing something).

> 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?

This would be one of the main points for running multiple simultaneous outcomes in one model, wouldn't it?

In my view, one would want to have variance-covariance Matrices (or possibly correlation matrices) for all random terms in the model (i.e. random effects and residuals). For that all random terms should end up in Matrices with their number of columns corresponding to the number of outcomes. Whereas this is the case for the random effects - e.g. ranef (mv1.lme) - the residuals are just assumed to be independt resulting in one vector - e.g. resid (mv1.lme).

This is also reflected in the degrees of freedom (were we inclined to rely on them, which, of course, we are not ;-) in that the dfs in such a model result from the total number of rows in the data set and not the number of observations per outcome.

The way to specify correlation structures of the errors in lme is via the corStruct methods. At least with the methods available in the package, it does not seem to be possible to model such a 'block'-structure of covariance (the residuals of the different outcomes co-varying with each other). Some of the pdMat classes might do the trick for random effects but I would not see how these could easily be applied to the residuals.

> ## 2) residual error correlated across variables at different 
> time point? (e.g., AR1)

The following might do the trick:
mv4.lme <- update (mv3.lme,
                   correlation= corAR1 (form= ~ 1 | subject/y1.ind))

(In the more general case with more than two outcomes, y1.ind would need to be replaced by a factor indicating the differen outcomes, i.e. y1.ind, y2.ind, ... "de-dummied", so to speak).

> ## 3) different AR1 coefficient for different variables?

No idea for that, either.
 
> ## 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)

This includes the random slopes which were not included above, correct?

I am currently thinking about implementing such models in flavour of Bugs and just use the simple but fast methods in lme/lmer as plausibility checks for the more complex models.

Best wishes, Lorenz
- 
Lorenz Gygax
Federal Veterinary Office
Centre for proper housing of ruminants and pigs
CH-8356 Ettenhausen / Switzerland


> ###############
> ## 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)))