Multivariate Mixed model in nlme

Question

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)

Answer 1

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)