# Multivariate Mixed model in nlme

I have two correlated response variables ($y_1,y_2$) explained by the same covariate set ($x_1,x_2$), each with mean = 0.

I have a random grouping factor ('group') which is heteroskedastic ($group \sim N(0,\sigma^2_{u1})$ for $y_1$ and $group \sim N(0,\sigma^2_{u2})$ for $y_2$).

I would like to model the correlation between the responses through the random error which is also heteroskedastic. $\epsilon_{1} \sim N (0,\sigma^2_{e1})$, $\epsilon_2 \sim N(0,\sigma^2_{e2}),$ and $Cov(\epsilon_{1},\epsilon_2)=\sigma_{e1,e2}$.

So, the equations for the $j^{\text{th}}$ observation in the $i^{\text{th}}$ group look like this:

$$y_{1ij} = \beta_{11} x_{1ij} + \beta_{12}x_{2ij} + group_{1i} + \epsilon_{1ij} \\ y_{2ij} = \beta_{21} x_{1ij} + \beta_{22}x_{2ij} + group_{2i} + \epsilon_{2ij} \\$$

I am using the following code to generate an example dataset. It is followed by the nlme code that I attempted (which gave me an error). I would appreciate any guidance on how to do this right.

Thanks!

library(MASS)
library(matrixcalc)
require(reshape2)
library(nlme)

ni <- c(rep(15,6),rep(20,6))
smpl <- sum(ni)
des <- factor(x=rep(x=1:length(ni),times=ni))
Z   <- model.matrix(~des-1,data=des)

mydata <- data.frame(des)
colnames(mydata) <- c('group')

means <- c(0,0,0,0)

# Covariance parameter values of Responses and Covariates
y1_sig <- 7
y2_sig <- 80
x1_sig <- 100
x2_sig <- 150

cov_fe <- matrix(c(
7*y1_sig     ,  -0.4*y1_sig*y2_sig,    0.75*y1_sig*x1_sig,   -0.65*y1_sig*x2_sig,
-0.4*y2_sig*y1_sig,       y2_sig*y2_sig,   -0.6 *y2_sig*x1_sig,    0.8*y2_sig*x2_sig,
0.75*x1_sig*y1_sig,  -0.6*x1_sig*y2_sig,         x1_sig*x1_sig,   -0.5*x1_sig*x2_sig,
-0.65*x2_sig*y1_sig,   0.8*x2_sig*y2_sig,   -0.5 *x2_sig*x1_sig,        x2_sig*x2_sig
),4,4)

set.seed(101)
fe <- mvrnorm( n = smpl, mu = means, Sigma = cov_fe, tol = 1e-6, empirical = FALSE)
mydata$X1 <- fe[,3] mydata$X2 <- fe[,4]

# random values
u1_sig <- 10
e1_sig <- 8
u2_sig <- 150
e2_sig <- 100

cov_ref <- diag(c(u1_sig*u1_sig ,u2_sig*u2_sig   ))
cov_rer <- diag(c(e1_sig*e1_sig , e2_sig*e2_sig  ))

ref <- mvrnorm(n = length(ni), mu = rep(0,2), Sigma = cov_ref, tol = 1e-6, empirical = FALSE)
reff <- Z %*% ref

rer <- mvrnorm( n = smpl,  mu = rep(0,2), Sigma = cov_rer,  tol = 1e-6, empirical = FALSE)

mydata$Y1 <- fe[,1] + reff[,1] + rer[,1] mydata$Y2 <- fe[,2] + reff[,2] + rer[,2]

# Stacking data into a univariate form
resp1 <- data.frame(mydata$Y1, mydata$Y2)
mydata_resp = melt(resp1,  variable_name='Y')

mat <- cbind(mydata$X1, mydata$X2, mydata$group) # New dataset mydata1 <- data.frame(direct.prod(diag(2),mat)) colnames(mydata1) <- c('X11','X21','group1','X12','X22','group2') mydata1$variable <- mydata_resp$variable mydata1$value <- mydata_resp$value mydata1$group1 <- factor(mydata1$group1) # Create binary variable to differentiate between stacked responses mydata1$D       <- as.integer(mydata1$variable == "mydata.Y1") # nlme call mdl.nlme <- lme(fixed = value ~ 0 + X11 + X21 + X12 + X22 , random = ~ -1| D /group1, correlation = corSymm(form=~1|D), data=mydata1)  ## 1 Answer I know this is a super old post. But I had the same question so I used your codes to explore a little bit. Not sure if this is what you wanted. mydata1$$X1 <- with(mydata1, X11+X12) mydata1$$X2 <- with(mydata1, X21+X22) mydata1$$group <- as.numeric(mydata1$$group1) + as.numeric(mydata1$group2)

mdl.nlme <- lme(fixed = value ~ -1 + D + X1 + X2 + X1:D + X2:D, random = ~ -1 + D | group, data=mydata1)
summary(mdl.nlme)