5
$\begingroup$

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)
$\endgroup$

1 Answer 1

1
$\begingroup$

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)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.