使用miceadds的ml.lmer方法插补缺失数据报错及交互项疑问

1. 数据与插补代码 #

待插补数据集 #

num = structure(list(f = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4), q = c(0, 0, 0, 1, 1, 1, 2, 2, 2, 0, 0, 0, 1, 1, 1, 2, 2, 2, 0), w = c(0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0), m = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), VAR1 = c(-4.95682177429289, -1.71885298325925, -1.87583453393747, -10.8948387329942, -8.08080625559951, -5.91060080742435, -0.441608813541267, -2.99216154761941, 2.84774625239502, NA, NA, NA, 9.99755147410283, 4.43052737120724, 14.5484885681107, 2.76708477891915, -1.16534086049291, 9.89837800524645, -0.532905309286036), VAR2 = c(1.97016982613686, 3.56723371979529, 3.51987609038002, -4.06622139517368, -2.12658955909922, -0.170511738732546, 7.61811332102057, 3.13322534451211, 9.57371767396202, NA, NA, NA, 11.984901005381, 5.72441513178024, 16.4228400177422, 6.70020469959361, 1.01261439823686, 11.6615475503401, 3.77633131099577), VAR3 = c(10.1560372560437, 10.1801786007975, 10.3578648604893, -0.140643528591531, 2.35612280120139, 5.07415319782397, 16.1107541935183, 9.04452647025277, 16.0511108882419, NA, NA, NA, 12.8200201578453, 6.83486757060971, 15.3398531136545, 8.2812335284254, 2.23463772852383, 13.3391933488311, 8.0877295387725), VAR4 = c(0.498180346940768, 3.00685583337094, 2.35322477767466, 0.282762095301972, 0.855921468843172, 2.66418197181061, 2.38290312756355, 2.67461144477413, 2.62873089735811, NA, NA, NA, 7.40263746337343, 5.65903049181366, 8.10425764744729, 3.64442942165172, 0.854005199532879, 4.05097634729794, 3.18750289721644), VAR5 = c(5.07634106445921, 5.4076172462957, 3.52593010658109, 3.29635357505705, 3.23171250638432, 5.13231721794768, 5.82362215060103, 4.50936077199033, 3.73864718826752, NA, NA, NA, 6.79332702092993, 4.99273093417307, 8.33859938657718, 5.99653923147157, 1.73552155550579, 4.99902596680605, 4.1363029065247), VAR6 = c(7.77422356168883, 8.82730417906284, 6.36866211542873, 4.85946360329115, 4.99276643284351, 7.11584125477072, 9.43260705283248, 6.17689370466653, 5.8380426743645, NA, NA, NA, 6.36662668308545, 4.75899543574918, 7.1874282896892, 6.15338369358482, 1.42804547517576, 5.7142481318893, 3.55850893758603)), row.names = c(NA, -19L), class = c("tbl_df", "tbl", "data.frame"))

原始插补代码 #

library(mice)
library(miceadds) 

predMatrix1 <- make.predictorMatrix(data = num)
impMethod1 <- make.method(data = num)
names = names(num[,5:10])

# method for lower-level variables (x, y, and z)
impMethod1[names] <- "ml.lmer"

predMatrix[, c("f", "q", "w", "m")] <- 0

level <- character(ncol(num))
names(level) <- colnames(num)

# 关注q和w的交互项，不清楚是否需要设置ID或特殊处理交互
vec = "w"
for (j in names) {
  level[j] <- vec
}

cluster <- list()
for (i in names) {
  cluster[[i]] <- c("q")
}

imputed = mice(num, method = impMethod1, predictorMatrix = predMatrix1, maxit = 1,
     m = 20, levels_id = cluster, variables_levels = level)

2. 报错信息 #

iter imp variable
  1   1  VAR1Error in if (stats::sd(weight.obs) > 0) { : 
  missing value where TRUE/FALSE needed 

3. 报错原因及修复方案 #

核心原因 #

1. 预测矩阵变量名错误：代码定义了predMatrix1，但后续修改时误用未定义的predMatrix，导致实际传入mice的预测矩阵未正确调整，逻辑混乱。
2. 关键变量被错误屏蔽：ml.lmer依赖聚类变量（q）和水平变量（w）构建混合效应模型，代码尝试将这些变量在预测矩阵中设为0，导致模型无法获取分组信息，计算观测权重标准差时出现缺失值，触发条件判断错误。
3. 变量类型不匹配：q和w为数值型，但作为分组/水平变量，应转换为因子型，否则混合模型无法正确识别其分类属性。

修复后的代码 #

library(mice)
library(miceadds) 

# 转换分类变量为因子
num$q <- as.factor(num$q)
num$w <- as.factor(num$w)

predMatrix1 <- make.predictorMatrix(data = num)
impMethod1 <- make.method(data = num)
var_names <- names(num[,5:10])

# 设置插补方法为ml.lmer
impMethod1[var_names] <- "ml.lmer"

# 仅屏蔽不需要的变量f和m，保留q和w作为预测变量
predMatrix1[, c("f", "m")] <- 0

# 定义每个插补变量的水平变量和聚类变量
variables_levels <- setNames(rep("w", length(var_names)), var_names)
levels_id <- setNames(rep(list("q"), length(var_names)), var_names)

# 执行插补（建议将maxit设为5以保证收敛）
imputed <- mice(num, method = impMethod1, predictorMatrix = predMatrix1, 
                maxit = 5, m = 20, levels_id = levels_id, 
                variables_levels = variables_levels)

4. 关于variables_levels处理交互项的说明 #

variables_levels仅用于指定嵌套结构中的水平变量，无法直接处理q和w的交互项。若需纳入交互项，需手动创建交互项列并加入预测矩阵：

# 创建q和w的交互项
num$q_w <- interaction(num$q, num$w, sep = "_")

# 更新预测矩阵，允许交互项作为预测变量
predMatrix1[, "q_w"] <- 1

更新后即可在插补模型中使用该交互项作为预测变量，提升模型拟合效果。

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