[R-sig-ME] Multivariate mixed effects model
Afshartous, David
DAfshartous at med.miami.edu
Mon Jul 6 17:59:32 CEST 2009
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
More information about the R-sig-mixed-models
mailing list