#StackBounty: #r #mixed-model #lme4-nlme Mixed effects model: Compare random variance component across levels of a grouping variable

Bounty: 100

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 slope variance in the groups defined by condition (i.e., the variance component highlighted in red, evaluated at condition=”control” vs. at condition=”experimental”). How would I do this (preferably in R)?

enter image description here


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 in 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:

  1. Is the question, “are the random intercepts 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?
  2. 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.


Get this bounty!!!

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.