线性混合模型模拟报错：Z矩阵与u矩阵维度不匹配

错误原因分析 #

出现non-conformable arguments错误的核心是Z矩阵维度与随机效应u的维度不匹配，具体问题如下：

• 原模型(1|Subject) + (1|Sex)表示：每个Subject水平对应一个随机截距，每个Sex水平对应一个随机截距。因此随机效应的总数是Subject的水平数（27个） + Sex的水平数（2个）=29个，而非代码中模拟的2个。
• 你构造的Z矩阵是model.matrix(~ -1 + Subject + Sex, Orthodont)，生成的是108行（观测数）×29列（27个Subject哑变量+2个Sex哑变量）的矩阵，但模拟的u只有2行（对应2个随机效应），两者相乘时维度不兼容。
• 对应的G矩阵也错误：原代码是2×2的对角矩阵，实际需要29×29的对角矩阵，其中前27个元素为sigma2_Subject，后2个为sigma2_Sex。

修正后的代码 #

1. 完整修正代码 #

library(lme4)
library(MASS)

# 拟合原模型并提取参数
data(Orthodont)
model_example <- lmer(distance~age + (1|Subject) + (1|Sex), Orthodont) 

# 规范提取方差成分与固定效应
vc <- VarCorr(model_example)
sigma2_Subject <- as.numeric(vc[["Subject"]])
sigma2_Sex <- as.numeric(vc[["Sex"]])
sigma_error <- sigma(model_example)^2
beta_int <- fixef(model_example)[1]
beta_age <- fixef(model_example)[2]

# 修正后的模拟函数
nboot = 100
sim <- function(nboot,
                dados,
                beta.int,
                beta.age,
                sigma2,
                sigma2_Subject, 
                sigma2_Sex  
){
  # 固定效应设计矩阵（基于输入数据集）
  X <- model.matrix(~ age, dados)
  
  # 固定效应系数向量
  beta <- c(beta.int, beta.age)
  
  # 构造随机效应的Z矩阵：拆分Subject和Sex的指示矩阵后合并
  Z_subj <- model.matrix(~ -1 + Subject, dados)
  Z_sex <- model.matrix(~ -1 + Sex, dados)
  Z <- cbind(Z_subj, Z_sex)
  
  # 随机效应协方差矩阵G：对应每个随机效应的方差
  n_subj <- ncol(Z_subj)
  n_sex <- ncol(Z_sex)
  G <- diag(c(rep(sigma2_Subject, n_subj), rep(sigma2_Sex, n_sex)))
  
  # 模拟误差项
  n_obs <- nrow(dados)
  e <- t(mvrnorm(n = nboot, 
                 mu = rep(0, n_obs),
                 Sigma = sigma2 * diag(n_obs),
                 empirical = FALSE))
  
  # 模拟随机效应：数量与Z矩阵列数一致
  u <- t(mvrnorm(n = nboot, 
                 mu = rep(0, n_subj + n_sex),
                 Sigma = G,
                 empirical = FALSE))
  
  # 计算模拟响应变量
  y <- X %*% beta %*% t(rep(1, nboot)) + Z %*% u + e
  return(y)
}

# 运行模拟
Y <- sim(nboot = 100,
         dados = Orthodont,
         beta.int = beta_int,
         beta.age = beta_age,
         sigma2 = sigma_error,
         sigma2_Subject = sigma2_Subject,
         sigma2_Sex = sigma2_Sex
)

# 组合模拟数据集
data_simulated <- cbind(Orthodont[,c("age", "Subject", "Sex")], Y)

2. 关键修正点说明 #

• Z矩阵构造：拆分Subject和Sex的指示矩阵后合并，确保每一列对应一个独立的随机效应（每个Subject水平、每个Sex水平各一列）。
• G矩阵调整：根据随机效应的总数生成对角矩阵，前27个元素对应Subject的随机截距方差，后2个对应Sex的随机截距方差。
• 随机效应u的模拟：模拟的随机效应数量与Z矩阵的列数完全匹配（27+2=29个），解决矩阵乘法的维度冲突。
• 代码健壮性优化：函数内直接使用输入的dados构造设计矩阵，避免依赖外部对象；用fixef()、sigma()等官方函数提取参数，比直接访问对象内部属性更稳定。

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