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 recover the two correlations `r1` and `r2` through variance partitioning or parameterization?**