How can we simulate correlated random variables that vary at different levels in a multilevel/mixed effects setting?

Question

I am very familiar with generating correlated random variables from a multivariate normal distribution.

This question is about doing that in a multilevel setting, where variables only vary at particular levels of a grouping variable.

So let's say we have one grouping factor - group and we wish to simulate 2 random variables, x1 and x2 where x1 varies at the lower level, and x2 varies at the upper level. Let's say the lower level has 100 (n1 = 100) observations, and the upper level has 20 distinct observations (n2 = 20, but obviously each value is replicated 5 times so that the group sizes are all equal (ie 5 per group).

How can we simulate x1 and x2 so that sd(x1) = 5 and sd(x2) = 3 and Cov(x1,x2) = 2 ?

I do not need any code. I have some code but would like some feedback about the approach, which is as follows:

In this method of simulating data for a two-level hierarchical model, we aim to generate two random variables, x1 and x2, where x1 varies at the lower (individual) level and x2 varies at the upper (group) level. The key goal is to ensure that the standard deviations of x1 and x2 are specified as 5 and 3 , respectively, and that the covariance between the expanded version of x2 (replicated across individuals within groups) and x1 is set to 2 (for the purposes of this simulation). The process begins by generating x2 at the group level with 20 distinct values for 20 groups, each with a standard deviation of 3. These group-level values are then replicated across individuals within each group to create the expanded x2. To achieve the desired covariance between x1 and x2, we calculate a shared component that correlates x1 with the expanded x2. This shared component is derived from the covariance formula, dividing the desired covariance by the variance of the expanded x2. We then generate x1 by adding independent noise to the shared component, ensuring that the overall variance of x1 equals 5. This method attempts to allow for controlled generation of correlated data across hierarchical levels, ensuring the correct standard deviations and covariance. This is my code that implements this.

# Parameters
n1 <- 1000  # Total number of observations at level 1 (individual level)
n2 <- 200   # Total number of groups at level 2 (group level)
group_size <- n1 / n2  # Size of each group

# Desired standard deviations and covariance
sd_x1 <- 5
sd_x2 <- 3
cov_x1_x2 <- 2

# Number of simulations
n_sim <- 100

vec_sd_1 <- numeric(n_sim)
vec_sd_2 <- numeric(n_sim)
vec_cov <- numeric(n_sim)

set.seed(15)

for (i in 1:n_sim) {

  # 1. Generate group-level variable x2 (with sd = 3)
 x2_group <- rnorm(n2, mean = 0, sd = sd_x2)

 # 2. Replicate x2 for each individual in the group (expand x2)
 x2 <- rep(x2_group, each = group_size)  # This makes x2 length n1

 # 3. Calculate the correct shared component based on expanded x2
 shared_component <- cov_x1_x2 / var(x2)  # The part of x1 that is correlated with the expanded x2

 # 4. Generate individual-level random variable x1
 x1 <- x2 * shared_component + rnorm(n1, mean = 0, sd = sqrt(sd_x1^2 - shared_component^2 * var(x2)))

 # Check results
 group <- rep(1:n2, each = group_size)

 # Data frame with simulated data
 data <- data.frame(group = factor(group), x1 = x1, x2 = x2)

 vec_sd_1[i]<- sd(data$x1)
 vec_sd_2[i] <- sd(data$x2)
 vec_cov[i] <- cov(data$x1, data$x2) 
 
}

mean(vec_sd_1) 
mean(vec_sd_2)
mean(vec_cov)

And these were the results:

> mean(vec_sd_1) 
[1] 4.990359
> mean(vec_sd_2)
[1] 2.994848
> mean(vec_cov)
[1] 2.003473

So it looks good, any feedback would be appreciated !

Answer 1

I like your approach, @Robert Long (+1), and my response should not be viewed as an answer. You are highlighting an important aspect of hierarchical data and the models we use to analyze them. A variable measured at the lower level (x1) is limited in how it can interact with a variable measured at the higher level (x2). The limitation comes from the fact that x1 has both within- and between-group variance, and if you want it to correlate with higher level variable x2, you have to do so only using the proportion of its total variance that sits at the higher level.

When I teach MLM and even in some of my responses on CV, I point out that level 1 predictors have within and between variation, which can be quantified by an ICC. This has implications for locating the level at which a predictor can have the most influence on the outcome. So in constructing your simulation, you had to account for this fact that x1 and x2 are correlated at the between level, not within. These lines of code do the trick elegantly:

# 3. Calculate the correct shared component based on expanded x2
shared_component <- cov_x1_x2 / var(x2)  # The part of x1 that is correlated with the expanded x2

# 4. Generate individual-level random variable x1
x1 <- x2 * shared_component + rnorm(n1, mean = 0, sd = sqrt(sd_x1^2 - shared_component^2 * var(x2)))

I am not sure I would do it differently, or that I could have coded it up as you have. I thought about possibly generating x1 at the between level by incorporating it's correlation w/ x2. Then one would have to add the noise at the within level, ensuring that the standard deviation of the noise at the within and between levels add up to 5. This is very similar to what you did, but again, yours is cleaner. Great work and thanks for posting your solution here for others to see!