I want to retrieve the correlations in a multivariate dataset. Let me first start with a simple case with three variables among which the first two are correlated. In other words, the three variables are assumed to follow a trivariate Gaussian $(y_{1}, y_{2}, y_{3})' \sim N((0, 0, 0)', S)$:

require(MASS)
r1 <- 0.5               # correlation value to be recovered
ns <- 2000              # number of samples

S  <- matrix(c(1,r1,0,  # correlation structure of trivariate data
               r1,1,0,
               0,0,1), nrow=3, ncol=3)

# simulated trivariate data
dat <- data.frame(f = c(rep(paste0('P',1:ns), 2), paste0('U',1:ns)),
                  y = c(mvrnorm(n=ns, mu=c(0, 0, 0), Sigma=S)))

In the data frame 'dat', each pair of samples from the two correlated variables (with a coefficient r1) are coded with the same label in the factor 'f'. Now we can construct the following model $y_{ij} = \alpha + \mu_i + \epsilon_{ij}$,

require(lme4)
m1 <- lmer(y ~ 1 + (1|f), data=dat)

With the variances from the model m1 output (each simulated dataset may lead to slightly different results):

summary(m1)

Random effects:
 Groups   Name        Variance Std.Dev.
 f        (Intercept) 0.4916   0.7012  
 Residual             0.5006   0.7075  

we can successfully recover the correlation r1 as below (which should be very close to the simulated r1 value, 0.5): $\frac{0.4916}{0.4916+0.5006} \approx 0.5$

tmp <- unlist(lapply(VarCorr(m1), `[`, 1))
# recover the correlation r1
tmp/(tmp+attr(VarCorr(m1), "sc")^2)

Now Let's switch to a case with 5 variables among which the first and second as well as the third and fourth are correlated. In other words, the five variables are assumed to follow a pentavariate Gaussian $(y_{1}, y_{2}, y_{3}, y_{4}, y_{5})' \sim N((0, 0, 0, 0, 0)', S)$:

r1 <- 0.2; r2 <- 0.8        # correlation value to be recovered
ns <- 2000                  # number of samples
S  <- matrix(c(1,r1,0,0,0,  # correlation structure of pentavariate data
               r1,1,0,0,0,  # the first and second variables are correlated
               0,0,1,r2,0,  # the third and fourth variables are correlated
               0,0,r2,1,0,
               0,0,0,0,1), nrow=5,ncol=5)

# simulated data
dat <- data.frame(f = c(rep(paste0('P',1:ns), 2), rep(paste0('T',1:ns), 2), paste0('S',1:ns)),
                  R=c(rep('P',2*ns), rep('T',2*ns), rep('U', ns)),
                  y = c(mvrnorm(n=ns, mu=rep(0,5), Sigma=S)))

In the data frame dat, the first and second variables are correlated (with coefficient r1); each pair of their samples are coded together with the same label in the factor f. Similarly, the third and fourth variables are correlated (with coefficient r2); each pair of their samples are coded together with the same label in the factor f. Because of the correlation structure, all the samples are categorized into 3 (instead of 5) levels in the factor R. Our goal is to use the simulated data to recover r1 and r2.

With the following dummy coding

dat$R1 <- as.numeric(dat$R=='P')   # dummy code the first and second variables (r1)
dat$R2 <- as.numeric(dat$R=='T')   # dummy code the third and fourth variables (r2)

I have considered the following two possible models

m2 <- lmer(y ~ 1 + (0+R1|f) + (0+R2|f), data=dat)
m3 <- lmer(y ~ 1 + (1|f) + (0+R2|f), data=dat)

To avoid over-parameterization, I don't include a term (1|f) (or (0+R1|f)) in the model m2 (or m3). However, I've been struggling to figure out a way to recover the correlations r1 and r2 with the variances from the random effects based on the models m2 and m3. In other words, the variances from the random effects in m2 and m3 don't seem to allow me to reconstruct r1 and r2.

summary(m2)

Random effects:
 Groups   Name Variance Std.Dev.
 f        R1   0.2731   0.5226  
 f.1      R2   0.5374   0.7331  
 Residual      0.6718   0.8196  

and

summary(m3)

Random effects:
 Groups   Name        Variance Std.Dev.
 f        (Intercept) 0.4004   0.6327  
 f.1      R2          0.2167   0.4655  
 Residual             0.5125   0.7159 

Of course we could use a workaround solution by reducing the situation into two cases with three variables:

# workaround solution
m4 <- lmer(y ~ 1 + (1|f), data=dat[dat$R2!=1,])	# recover r1 like model m1 above
m5 <- lmer(y ~ 1 + (1|f), data=dat[dat$R1!=1,])    # recover r2 like model m1 above

Then we can adopt the same strategy as the first example with three variables and recover r1 and r2 separately. However, I would really want to find a way to recover r1 and r2 directly using the full data with models like m2 and m3.

The models m2 and m3 are not the right formulation for the retrieval of r1 and r2? So, I'm stuck: Is there a way to properly parameterize the effects so that the two correlations could be recovered?