R语言mice包二进制变量delta敏感性分析代码异常问询

背景 #

主分析前使用mice包执行多重插补，结局变量为功能衰退（FunctionalDecline），预测变量包含连续/二元变量：连续变量采用pmm插补，二元变量采用logreg插补。针对结局变量的缺失值开展delta调整敏感性分析，假设缺失值更易对应衰退者，取delta=log(4.3)（参考Rezvan等人研究），基于Leacy等人论文自定义了mice.impute.logreg.sens插补函数。

自定义插补函数代码 #

library(mice)
mice.impute.logreg.sens <- function (y, ry, x, delta, ...) {
  x <- cbind(1, as.matrix(x))
  expr <- expression(glm.fit(x[ry, ], y[ry], family = binomial(link = logit)))
  fit <- suppressWarnings(eval(expr))
  fit.sum <- summary.glm(fit)
  beta <- coef(fit)
  beta[1] <- beta[1] + delta
  rv <- t(chol(fit.sum$cov.unscaled))
  beta.star <- beta + rv %*% rnorm(ncol(rv))
  p <- 1 / (1 + exp(-(x[!ry, ] %*% beta.star)))
  vec <- (runif(nrow(p)) <= p)
  vec[vec] <- 1
  if (is.factor(y)) {
    vec <- factor(vec, c(0, 1), levels(y))
  }
  return (vec)
}

数据集应用代码 #

myvars <- c("FunctionalDecline", "Age", "Sex", "GDS4")
realdataset <- dataset[myvars]
realdataset$FunctionalDecline <- factor(realdataset$FunctionalDecline)
realdataset$Sex <- factor(realdataset$Sex)
ini <- mice(realdataset, maxit = 0)
meth <- ini$meth
meth["FunctionalDecline"] = "logreg.sens"
imputeddataset <- mice(realdataset, meth = meth, seed = 3, m = 5, maxit = 10, delta = log(4.3), print = F)
summary(pool(with(imputeddataset, glm(FunctionalDecline~Age + Sex + GDS4, family = binomial))))

异常现象 #

• 当delta=log(1.0)（无调整）时，logistic回归结果与未做delta调整的插补数据集回归结果不一致；
• delta=log(4.3)调整后，所有预测变量均失去显著性，但“粗略最坏情景”（将缺失的功能衰退值全设为1）分析中多数预测变量仍显著；
• 当delta取log(280)（对应缺失值中99%为衰退者）时，结果与上述“粗略最坏情景”完全不符。
  调整迭代次数、插补数量或访问序列后，结果无实质变化。

错误原因分析 #

1. delta参数传递失效 #

mice函数的额外参数（如delta）不会自动传递给自定义插补函数，当前代码中mice.impute.logreg.sens的delta参数无法正确获取mice()中传入的值，导致插补时delta未按预期生效，直接引发delta=log(1.0)时结果与标准插补不一致的问题。

2. 截距调整逻辑偏差 #

Leacy等人的方法中，delta对应**优势比（OR）**的调整，映射到logit尺度是对截距的调整，但当前代码存在两个问题：

• 未考虑因子变量的水平顺序：若FunctionalDecline的参考水平不是0，直接调整截距会导致概率预测方向错误；
• 仅单纯叠加delta到原始截距，未结合敏感性分析的假设（缺失值更易为衰退者）确认调整方向的正确性。

3. 协方差矩阵使用错误 #

代码中fit.sum$cov.unscaled是原始模型的协方差矩阵，调整截距后未更新协方差结构，导致生成的beta.star（带随机误差的系数）不符合敏感性分析的假设分布，进而使得插补概率偏离预期，极端delta值下无法趋近于“全设为1”的情景。

4. 因子变量处理漏洞 #

当y是因子时，vec <- factor(vec, c(0, 1), levels(y))的逻辑存在问题：如果levels(y)的顺序不是c(0,1)，会导致插补值的因子水平映射错误，插补结果的实际含义与预期不符。

修正建议 #

1. 正确传递delta参数：在mice()调用中通过...传递delta，并在自定义函数内通过args <- list(...)提取，例如：

   mice.impute.logreg.sens <- function (y, ry, x, ...) {
     args <- list(...)
     delta <- args$delta
     # 后续逻辑
   }
   

2. 修正截距调整逻辑：先确认y的水平顺序，确保delta的调整方向符合“缺失值更易为衰退者”的假设，例如若衰退对应y=1，则增加截距以提高预测为1的概率；
3. 更新协方差结构：调整截距后，基于新的系数重新计算协方差矩阵，或基于mice内置logreg函数的逻辑修改截距部分；
4. 规范因子变量处理：根据y的实际水平映射插补值，例如：

   if (is.factor(y)) {
     vec <- factor(ifelse(vec, levels(y)[which(levels(y)=="1")], levels(y)[which(levels(y)=="0")]), 
                   levels = levels(y))
   }
   

5. 验证插补结果：插补后检查缺失值的插补比例，例如当delta=log(280)时，插补为1的比例应接近99%，以此验证函数是否正确实现假设。

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