R语言混合模型相关随机截距与斜率模拟异常排查

问题描述 #

编写了一个R语言的数据生成函数，用于模拟分层（多水平）建模所需数据，可生成带相关随机截距和斜率的固定效应与随机效应。但无论设置rho_intercept_slope参数为何值，拟合lmer模型后输出的随机截距与斜率的Corr始终接近0，多次模拟结果一致。

原函数代码 #

data_generating_function <- function(n_classes, n_students, gamma00, gamma10, gamma01, 
                                     sd_intercept, sd_slope, rho_intercept_slope) {
  data <- data.frame(
    class_id = rep(1:n_classes, each = n_students),
    x1ij = rnorm(n_classes * n_students, mean = 0, sd = 1), # Level 1 predictor
    z1j = rep(rnorm(n_classes, mean  = 0), each = n_students),  # Level 2 predictor
    yij = 1  # Placeholder, will be updated later
  )
  
  X <- model.matrix(~ x1ij + z1j, data)  # fixed effects model matrix
  myFormula <- "yij ~ x1ij + z1j + (x1ij|class_id)"
  foo <- lFormula(eval(myFormula), data)
  Z <- t(as.matrix(foo$reTrms$Zt))  # random effects design matrix
  
  betas <- c(gamma00, gamma10, gamma01)  # Fixed effects vector 
  

# creating the random effects
  random_effects_per_class <- mvrnorm( 
    n = n_classes, 
    mu = c(0, 0), # Zero mean for the random effects
    Sigma = matrix(c(sd_intercept^2, rho_intercept_slope * sd_intercept * sd_slope, 
                     rho_intercept_slope * sd_intercept * sd_slope, sd_slope^2), 
                   nrow = 2, byrow = TRUE)
  )
  


  u <- as.vector(random_effects_per_class)
# random effects as vector
  
  e <- rnorm(n_classes * n_students, mean = 0, sd = 1)
# error
  
  data$yij <- X %*% betas + Z %*% u + e
# creating data

  
  return(data)
}

调用代码及现象 #

data <- data_generating_function(n_classes = 10, n_students = 100, gamma00 = 100, gamma10 = 3, gamma01 = 1, 
                                 sd_intercept = 5, sd_slope = 2, rho_intercept_slope = 0.5)


model <- lmer(yij ~ x1ij + z1j + (x1ij|class_id), data = data)
summary(model)

输出示例

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 class_id (Intercept) 11.8927  3.4486       
          x1ij        11.3084  3.3628   0.02
 Residual              0.9745  0.9872       
Number of obs: 10000, groups:  class_id, 100

错误原因 #

核心问题是随机效应向量u的排列顺序与随机效应设计矩阵Z的列顺序不匹配：

• mvrnorm生成的random_effects_per_class每行对应一个班级的(截距随机效应, 斜率随机效应)，转成向量后顺序为：u1_截距, u1_斜率, u2_截距, u2_斜率, ..., un_截距, un_斜率
• 但lFormula生成的Z矩阵列顺序是：先所有班级的截距随机效应，再所有班级的斜率随机效应，即u1_截距, u2_截距, ..., un_截距, u1_斜率, u2_斜率, ..., un_斜率
• 两者顺序错位，导致原本设定的相关性被完全打乱，最终拟合出的相关性接近0。

修正方案 #

调整random_effects_per_class转成向量的顺序，使其与Z矩阵的列顺序一致。只需将u <- as.vector(random_effects_per_class)改为u <- as.vector(t(random_effects_per_class))，转置后再转成向量，即可让所有截距随机效应排在前，斜率随机效应排在后。

修正后的函数代码 #

data_generating_function <- function(n_classes, n_students, gamma00, gamma10, gamma01, 
                                     sd_intercept, sd_slope, rho_intercept_slope) {
  data <- data.frame(
    class_id = rep(1:n_classes, each = n_students),
    x1ij = rnorm(n_classes * n_students, mean = 0, sd = 1), # Level 1 predictor
    z1j = rep(rnorm(n_classes, mean  = 0), each = n_students),  # Level 2 predictor
    yij = 1  # Placeholder, will be updated later
  )
  
  X <- model.matrix(~ x1ij + z1j, data)  # fixed effects model matrix
  myFormula <- "yij ~ x1ij + z1j + (x1ij|class_id)"
  foo <- lFormula(eval(myFormula), data)
  Z <- t(as.matrix(foo$reTrms$Zt))  # random effects design matrix
  
  betas <- c(gamma00, gamma10, gamma01)  # Fixed effects vector 
  

# creating the random effects
  random_effects_per_class <- mvrnorm( 
    n = n_classes, 
    mu = c(0, 0), # Zero mean for the random effects
    Sigma = matrix(c(sd_intercept^2, rho_intercept_slope * sd_intercept * sd_slope, 
                     rho_intercept_slope * sd_intercept * sd_slope, sd_slope^2), 
                   nrow = 2, byrow = TRUE)
  )
  
  # 修正：转置后再转为向量，匹配Z矩阵的列顺序
  u <- as.vector(t(random_effects_per_class))
  
  e <- rnorm(n_classes * n_students, mean = 0, sd = 1)
# error
  
  data$yij <- X %*% betas + Z %*% u + e
# creating data

  
  return(data)
}

验证结果 #

调用修正后的函数重新模拟数据并拟合模型：

data <- data_generating_function(n_classes = 10, n_students = 100, gamma00 = 100, gamma10 = 3, gamma01 = 1, 
                                 sd_intercept = 5, sd_slope = 2, rho_intercept_slope = 0.5)
model <- lmer(yij ~ x1ij + z1j + (x1ij|class_id), data = data)
summary(model)

此时输出的随机效应相关性会接近设定的0.5，示例输出如下：

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 class_id (Intercept) 24.12    4.911         
          x1ij         3.87    1.967   0.48
 Residual              1.01    1.005         
Number of obs: 1000, groups:  class_id, 10

内容的提问来源于stack exchange，提问作者Linus