Multivariate mixed model using nlme

Question

I have two correlated response variables (y1,y2) explained by the same covariate set (x1,x2), each with mean = 0.

I have a random grouping factor ('group') which is heteroskedastic (group ~ N(0,sigma_u1) for y1 and N(0,sigma_u2) for y2).

I would like to model the correlation between the responses through the random error which is also heteroskedastic.

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.

I apologize in advance if my code seems clunky.

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|group1/D,
                correlation = corSymm(form=~1|D),
                data=mydata1)