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 correlation value of r1 in the variance-covariance matrix S) are coded with the same label P in the factor 'f' while the samples for the third variable are coded with the label U. 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 4 variables among which the first and second as well as the third and fourth are correlated. In other words, the four variables are assumed to follow a quadrivariate Gaussian $(y_{1}, y_{2}, y_{3}, y_{4})' \sim N((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,    # correlation structure of quadrivariate data
               r1,1,0,0,    # the first and second variables are correlated
               0,0,1,r2,    # the third and fourth variables are correlated
               0,0,r2,1), nrow=4,ncol=4)

# simulated data
dat <- data.frame(f = c(rep(paste0('P',1:ns), 2), rep(paste0('T',1:ns), 2)),
                  R=c(rep('P',2*ns), rep('T',2*ns)),
                  y = c(mvrnorm(n=ns, mu=rep(0,4), 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 2 (instead of ) 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.3835   0.6193  
 f.1      R2   0.6537   0.8085  
 Residual      0.5012   0.7079

and

 summary(m3)

 Random effects:    
  Groups   Name        Variance Std.Dev.
  f        (Intercept) 0.3835   0.6193  
  f.1      R2          0.2702   0.5198  
  Residual             0.5012   0.7079 

Of course we could use a workaround solution by reducing the situation into two cases with two 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.

One noticeable aspect is that in the model m1 the total variance is conserved in the sense that the sum of the two variances from the modeling result ($0.4916+0.5006\approx 1$) is equal to the the variance (which is 1) of the simulating distribution $N((0, 0, 0)', S)$. In contrast, this is not true for the models m2 and m3: the 3 variances in summary(m2) and summary(m3) do not add up to 1. This may indicate that there is variance leakage happening in the two latter models? Or the models m2 and m3 are not parameterized properly for retrieving r1 and r2?

So, I'm stuck: Is there a way to properly parameterize the effects so that the two correlations could be recovered?