[R-sig-ME] Multivariate mixed effects model

Afshartous, David DAfshartous at med.miami.edu
Mon Jul 20 21:48:04 CEST 2009

Dear Lorenz and others interested,
Please see reply below.   In summary, the basic issue that remains is how to specify a general covariance matrix for the residual errors of several response variables at the same time point.
Any ideas much appreciated.

On 7/13/09 3:42 AM, "lorenz.gygax at art.admin.ch" <lorenz.gygax at art.admin.ch> wrote:

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.

<DA> The fact that the data is multivariate is actually not that much of an extension.  For instance, if we had only one response variable but say k groups with random effect variance or residual error variance stratified per group, the situation would be identical.  Just replace "response variable" by "group": Instead of k different response variables we have k groups for one variable and the data structure and model statement with appropriate dummy variables would be the same in lme/lmer.

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

<DA> one example would be to have the residual error of the different variables to be correlated at each time point as mentioned below.

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

<DA> Also, presumably the association structure of the random effects would often be of interest.   VarIdent allows different residual variances per variable, but does not allow them to be correlated.  Based on the list of varFunc classes (p.208, Pinheiro & Bates) I don't see how this would be done.  One wants something akin to pdSymm that is used for random effects.

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

<DA> Thanks.  I also noticed that the statement below accomplishes the same thing:
mv4.lme <- update (mv3.lme,   correlation= corAR1 ())

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

<DA> Yes

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

More information about the R-sig-mixed-models mailing list