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?**