Two residual covariance matrices for two blocks of observatiopns

Interfaces
1

I am trying to model y[1:(n1+n2),1:p] as independent multivariate normal where for the first set of 1:n1 observations, y has one covariance matrix and for the next set of (n1+1):(n1+n2), y has a different covariance matrix.

I tried to do this in brms using two sets of random intercepts. However, it does not appear to work very well, in terms of estimating the data generating Sigma_1 and Sigma_2. I understand that there is some non-identifiability for the diagonal variance terms in the analysis model.

Is there an easy way to specify this in brms?

# Data generation
y1 = mvrnorm(n1,rep(0,p),sigma_1)
y2 = mvrnorm(n2,rep(0,p),sigma_2)
y = rbind(y1,y2)

# Analysis code
ind_1 = c(rep(1,n1), rep(0,n2))
group_1= 1:(n1+n2) 
ind_2 = c(rep(0,n1), rep(1,n2))
group_2= 1:(n1+n2)

data_y = data.frame(y=y, ind_1=ind, ind_2=ind_2, 
                    group_1=group_1, group_2=group_2)

form_lst =list()
for (j in 1:p){
  form_lst[[j]] = bf(as.formula(paste("y.",j,"~(ind_1-1|c1|group_1) + (ind_2-1|c2|group_2)",sep="")))
  }

fit = brm(form_lst[[1]] + form_lst[[2]] + form_lst[[3]] + 
      form_lst[[4]] + form_lst[[5]] + form_lst[[6]] +  form_lst[[7]]+
     set_rescor(FALSE),
    data = data_y, chains = 4, cores = 4,
    iter = 20000, warmup = 5000)

Thank you.

2

I don’t know anything about brms, but you can easily do this in Stan itself with array slicing.