Mixed effects model: Compare random variance component across levels of a grouping variable Suppose I have $N$ participants, each of whom gives a response $Y$ 20 times, 10 in one condition and 10 in another.  I fit a linear mixed effects model comparing $Y$ in each condition.  Here's a reproducible example simulating this situation using the lme4 package in R:
library(lme4)
fml <- "~ condition + (condition | participant_id)"
d <- expand.grid(participant_id=1:40, trial_num=1:10)
d <- rbind(cbind(d, condition="control"), cbind(d, condition="experimental"))

set.seed(23432)
d <- cbind(d, simulate(formula(fml), 
                       newparams=list(beta=c(0, .5), 
                                      theta=c(.5, 0, 0), 
                                      sigma=1), 
                       family=gaussian, 
                       newdata=d))

m <- lmer(paste("sim_1 ", fml), data=d)
summary(m)

The model m yields two fixed effects (an intercept and slope for condition), and three random effects (a by-participant random intercept, a by-participant random slope for condition, and an intercept-slope correlation).
I would like to statistically compare the size of the by-participant random intercept variance across the groups defined by condition (i.e., compute the variance component highlighted in red separately within the control and experimental conditions, then test whether the difference in the size of the components is non-zero).  How would I do this (preferably in R)?


BONUS
Let's say the model is slightly more complicated: The participants each experience 10 stimuli 20 times each, 10 in one condition and 10 in another.  Thus, there are two sets of crossed random effects: random effects for participant and random effects for stimulus.  Here's a reproducible example:
library(lme4)
fml <- "~ condition + (condition | participant_id) + (condition | stimulus_id)"
d <- expand.grid(participant_id=1:40, stimulus_id=1:10, trial_num=1:10)
d <- rbind(cbind(d, condition="control"), cbind(d, condition="experimental"))

set.seed(23432)
d <- cbind(d, simulate(formula(fml), 
                       newparams=list(beta=c(0, .5), 
                                      theta=c(.5, 0, 0, .5, 0, 0), 
                                      sigma=1), 
                       family=gaussian, 
                       newdata=d))

m <- lmer(paste("sim_1 ", fml), data=d)
summary(m)

I would like to statistically compare the magnitude of the random by-participant intercept variance across the groups defined by condition.  How would I do that, and is the process any different from the one in the situation described above?

EDIT
To be a bit more specific about what I'm looking for, I want to know:


*

*Is the question, "are the conditional mean responses within each condition (i.e., random intercept values in each condition) substantially different from each other, beyond what we would expect due to sampling error" a well-defined question (i.e., is this question even theoretically answerable)?  If not, why not?

*If the answer to question (1) is yes, how would I answer it? I would prefer an R implementation, but I'm not tied to the lme4 package -- for example, it seems as though the OpenMx package has the capability to accommodate multi-group and multi-level analyses (https://openmx.ssri.psu.edu/openmx-features), and this seems like the sort of question that ought to be answerable in an SEM framework.

 A: One relatively straight-forward way could be to use likelihood-ratio tests via anova as described in the lme4 FAQ.
We start with a full model in which the variances are unconstrained (i.e., two different variances are allowed) and then fit one constrained model in which the two variances are assumed to be equal. We simply compare them with anova() (note that I set REML = FALSE although REML = TRUE with anova(..., refit = FALSE) is completely feasible).
m_full <- lmer(sim_1 ~ condition + (condition | participant_id), data=d, REML = FALSE)
summary(m_full)$varcor
 # Groups         Name                  Std.Dev. Corr  
 # participant_id (Intercept)           0.48741        
 #                conditionexperimental 0.26468  -0.419
 # Residual                             1.02677     

m_red <- lmer(sim_1 ~ condition + (1 | participant_id), data=d, REML = FALSE)
summary(m_red)$varcor
 # Groups         Name        Std.Dev.
 # participant_id (Intercept) 0.44734 
 # Residual                   1.03571 

anova(m_full, m_red)
# Data: d
# Models:
# m_red: sim_1 ~ condition + (1 | participant_id)
# m_full: sim_1 ~ condition + (condition | participant_id)
#        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
# m_red   4 2396.6 2415.3 -1194.3   2388.6                         
# m_full  6 2398.7 2426.8 -1193.3   2386.7 1.9037      2      0.386

However, this test is likely conservative. For example, the FAQ says:

Keep in mind that LRT-based null hypothesis tests are conservative
  when the null value (such as σ2=0) is on the boundary of the feasible
  space; in the simplest case (single random effect variance), the
  p-value is approximately twice as large as it should be (Pinheiro and
  Bates 2000).

There are several alternatives:


*

*Create an appropriate test distribution, which usually consists of a mixture of $\chi^2$ distributions. See e.g.,
Self, S. G., & Liang, K.-Y. (1987). Asymptotic Properties of Maximum Likelihood Estimators and Likelihood Ratio Tests Under Nonstandard Conditions. Journal of the American Statistical Association, 82(398), 605. https://doi.org/10.2307/2289471
However, this is quite complicated.

*Simulate the correct distribution using RLRsim (as also described in the FAQ).
I will demonstrate the second option in the following:
library("RLRsim")
## reparametrize model so we can get one parameter that we want to be zero:
afex::set_sum_contrasts() ## warning, changes contrasts globally
d <- cbind(d, difference = model.matrix(~condition, d)[,"condition1"])

m_full2 <- lmer(sim_1 ~ condition + (difference | participant_id), data=d, REML = FALSE)
all.equal(deviance(m_full), deviance(m_full2))  ## both full models are identical

## however, we need the full model without correlation!
m_full2b <- lmer(sim_1 ~ condition + (1| participant_id) + 
                   (0 + difference | participant_id), data=d, REML = FALSE)
summary(m_full2b)$varcor
 # Groups           Name        Std.Dev.
 # participant_id   (Intercept) 0.44837 
 # participant_id.1 difference  0.13234 
 # Residual                     1.02677 

## model that only has random effect to be tested
m_red <- update(m_full2b,  . ~ . - (1 | participant_id), data=d, REML = FALSE)
summary(m_red)$varcor
 # Groups         Name       Std.Dev.
 # participant_id difference 0.083262
 # Residual                  1.125116

## Null model 
m_null <- update(m_full2b,  . ~ . - (0 + difference | participant_id), data=d, REML = FALSE)
summary(m_null)$varcor
 # Groups         Name        Std.Dev.
 # participant_id (Intercept) 0.44734 
 # Residual                   1.03571 

exactRLRT(m_red, m_full2b, m_null)
# Using restricted likelihood evaluated at ML estimators.
# Refit with method="REML" for exact results.
# 
#   simulated finite sample distribution of RLRT.
#   
#   (p-value based on 10000 simulated values)
# 
# data:  
# RLRT = 1.9698, p-value = 0.0719

As we can see, the output suggests that with REML = TRUE we would have gotten exact results. But this is left as an exercise to the reader.
Regarding the bonus, I am not sure if RLRsim allows simultaneous testing of multiple components, but if so, this can be done in the same way. 

Response to comment: 

So it is true, then, that in general the random slope $\theta_X$ allows the random intercept $\theta_0$ to vary across levels of $X$?

I am not sure this question can receive a reasonable answer.  


*

*A random intercept allows an idiosyncratic difference in the overall level for each level of the grouping factor. For example, if the dependent variable is response time, some participants are faster and some are slower. 

*A random slope allows each level of the grouping factor an idiosyncratic effect of the factor for which random slopes are estimated. For example, if the factor is congruency, then some participants can have a higher congruency effect than others. 


So do random-slopes affect the random-intercept? In some sense this might make sense, as they allow each level of the grouping factor a completely idiosyncratic effect for each condition. In the end, we estimate two idiosyncratic parameters for two conditions. However, I think the distinction between the overall level captured by the intercept and the condition specific effect captured by the random slope is a important and then the random slope cannot really affect the random intercept. However, it still allows each level of the grouping factor an idiosyncratic separately for each level of the condition. 
Nevertheless, my test still does what the original question wants. It tests whether the difference in variances between the two conditions is zero. If it is zero, then the variances in both conditions are equal. In other words, only if there is no need for a random-slope is the variance in both conditions identical. I hope that makes sense.
A: Your model 
m = lmer(sim_1 ~ condition + (condition | participant_id), data=d)

already allows the across-subject variance in the control condition to differ from the across-subject variance in the experimental condition. This can be made more explicit by an equivalent re-parametrization:
m = lmer(sim_1 ~ 0 + condition + (0 + condition | participant_id), data=d)

The random covariance matrix now has a simpler interpretation:
Random effects:
 Groups         Name                  Variance Std.Dev. Corr
 participant_id conditioncontrol      0.2464   0.4963       
                conditionexperimental 0.2074   0.4554   0.83

Here the two variances are precisely the two variances you are interested in: the [across-subjects] variance of conditional mean responses in the control condition and the same in the experimental condition. In your simulated dataset, they are 0.25 and 0.21. The difference is given by
delta = as.data.frame(VarCorr(m))[1,4] - as.data.frame(VarCorr(m))[2,4]

and is equal to 0.039. You want to test if it is significantly different from zero.
EDIT: I realized that the permutation test that I describe below is incorrect; it won't work as intended if the means in control/experimental condition are not the same (because then the observations are not exchangeable under the null). It might be a better idea to bootstrap subjects (or subjects/items in the Bonus case) and obtain the confidence interval for delta.
I will try to fix the code below to do that.

Original permutation-based suggestion (wrong)
I often find that one can save oneself a lot of trouble by doing a permutation test. Indeed, in this case it is very easy to set up. Let's permute control/experimental conditions for each subject separately; then any difference in variances should be eliminated. Repeating this many times will yield the null distribution for the differences.
(I do not program in R; everybody please feel free to re-write the following in a better R style.)
set.seed(42)
nrep = 100
v = matrix(nrow=nrep, ncol=1)
for (i in 1:nrep)
{
   dp = d
   for (s in unique(d$participant_id)){             
     if (rbinom(1,1,.5)==1){
       dp[p$participant_id==s & d$condition=='control',]$condition = 'experimental'
       dp[p$participant_id==s & d$condition=='experimental',]$condition = 'control'
     }
   }
  m <- lmer(sim_1 ~ 0 + condition + (0 + condition | participant_id), data=dp)
  v[i,] = as.data.frame(VarCorr(m))[1,4] - as.data.frame(VarCorr(m))[2,4]
}
pvalue = sum(abs(v) >= abs(delta)) / nrep

Running this yields the p-value $p=0.7$. One can increase nrep to 1000 or so.
Exactly the same logic can be applied in your Bonus case.
A: There's more than one way to test this hypothesis. For example, the procedure outlined by @amoeba should work. But it seems to me that the simplest, most expedient way to test it is using a good old likelihood ratio test comparing two nested models. The only potentially tricky part of this approach is in knowing how to set up the pair of models so that dropping out a single parameter will cleanly test the desired hypothesis of unequal variances. Below I explain how to do that.
Short answer
Switch to contrast (sum to zero) coding for your independent variable and then do a likelihood ratio test comparing your full model to a model that forces the correlation between random slopes and random intercepts to be 0:
# switch to numeric (not factor) contrast codes
d$contrast <- 2*(d$condition == 'experimental') - 1

# reduced model without correlation parameter
mod1 <- lmer(sim_1 ~ contrast + (contrast || participant_id), data=d)

# full model with correlation parameter
mod2 <- lmer(sim_1 ~ contrast + (contrast | participant_id), data=d)

# likelihood ratio test
anova(mod1, mod2)

Visual explanation / intuition
In order for this answer to make sense, you need to have an intuitive understanding of what different values of the correlation parameter imply for the observed data. Consider the (randomly varying) subject-specific regression lines. Basically, the correlation parameter controls whether the participant regression lines "fan out to the right" (positive correlation) or "fan out to the left" (negative correlation) relative to the point $X=0$, where X is your contrast-coded independent variable. Either of these imply unequal variance in participants' conditional mean responses. This is illustrated below:

In this plot, we ignore the multiple observations that we have for each subject in each condition and instead just plot each subject's two random means, with a line connecting them, representing that subject's random slope. (This is made up data from 10 hypothetical subjects, not the data posted in the OP.)
In the column on the left, where there's a strong negative slope-intercept correlation, the regression lines fan out to the left relative to the point $X=0$. As you can see clearly in the figure, this leads to a greater variance in the subjects' random means in condition $X=-1$ than in condition $X=1$.
The column on the right shows the reverse, mirror image of this pattern. In this case there is greater variance in the subjects' random means in condition $X=1$ than in condition $X=-1$.
The column in the middle shows what happens when the random slopes and random intercepts are uncorrelated. This means that the regression lines fan out to the left exactly as much as they fan out to the right, relative to the point $X=0$. This implies that the variances of the subjects' means in the two conditions are equal.
It's crucial here that we've used a sum-to-zero contrast coding scheme, not dummy codes (that is, not setting the groups at $X=0$ vs. $X=1$). It is only under the contrast coding scheme that we have this relationship wherein the variances are equal if and only if the slope-intercept correlation is 0. The figure below tries to build that intuition:

What this figure shows is the same exact dataset in both columns, but with the independent variable coded two different ways. In the column on the left we use contrast codes -- this is exactly the situation from the first figure. In the column on the right we use dummy codes. This alters the meaning of the intercepts -- now the intercepts represent the subjects' predicted responses in the control group. The bottom panel shows the consequence of this change, namely, that the slope-intercept correlation is no longer anywhere close to 0, even though the data are the same in a deep sense and the conditional variances are equal in both cases. If this still doesn't seem to make much sense, studying this previous answer of mine where I talk more about this phenomenon may help.
Proof
Let $y_{ijk}$ be the $j$th response of the $i$th subject under condition $k$. (We have only two conditions here, so $k$ is just either 1 or 2.) Then the mixed model can be written
$$
y_{ijk} = \alpha_i + \beta_ix_k + e_{ijk},
$$
where $\alpha_i$ are the subjects' random intercepts and have variance $\sigma^2_\alpha$, $\beta_i$ are the subjects' random slope and have variance $\sigma^2_\beta$, $e_{ijk}$ is the observation-level error term, and $\text{cov}(\alpha_i, \beta_i)=\sigma_{\alpha\beta}$.
We wish to show that
$$
\text{var}(\alpha_i + \beta_ix_1) = \text{var}(\alpha_i + \beta_ix_2) \Leftrightarrow \sigma_{\alpha\beta}=0.
$$
Beginning with the left hand side of this implication, we have
$$
\begin{aligned}
\text{var}(\alpha_i + \beta_ix_1) &= \text{var}(\alpha_i + \beta_ix_2) \\
\sigma^2_\alpha + x^2_1\sigma^2_\beta + 2x_1\sigma_{\alpha\beta} &= \sigma^2_\alpha + x^2_2\sigma^2_\beta + 2x_2\sigma_{\alpha\beta} \\
\sigma^2_\beta(x_1^2 - x_2^2) + 2\sigma_{\alpha\beta}(x_1 - x_2) &= 0.
\end{aligned}
$$
Sum-to-zero contrast codes imply that $x_1 + x_2 = 0$ and $x_1^2 = x_2^2 = x^2$. Then we can further reduce the last line of the above to
$$
\begin{aligned}
\sigma^2_\beta(x^2 - x^2) + 2\sigma_{\alpha\beta}(x_1 + x_1) &= 0 \\
\sigma_{\alpha\beta} &= 0,
\end{aligned}
$$
which is what we wanted to prove. (To establish the other direction of the implication, we can just follow these same steps in reverse.)
To reiterate, this shows that if the independent variable is contrast (sum to zero) coded, then the variances of the subjects' random means in each condition are equal if and only if the correlation between random slopes and random intercepts is 0. The key take-away point from all this is that testing the null hypothesis that $\sigma_{\alpha\beta} = 0$ will test the null hypothesis of equal variances described by the OP.
This does NOT work if the independent variable is, say, dummy coded. Specifically, if we plug the values $x_1=0$ and $x_2=1$ into the equations above, we find that
$$
\text{var}(\alpha_i) = \text{var}(\alpha_i + \beta_i) \Leftrightarrow \sigma_{\alpha\beta} = -\frac{\sigma^2_\beta}{2}.
$$
