R mice package : error in 2l. norm 只有1个NA

Question

在调用 mice 时指定“2l.norm”方法后，我偶然发现了一条错误消息，变量仅包含 1 个 NA。考虑到这些变量的缺失数据量非常小，我意识到这是一个非常小的问题。但是，最好也考虑这些的数据结构。

我使用所有人都可以访问的数据库 ChickWeight 数据集重新创建了这种情况。我非常清楚这个问题也可能是我在执行程序时出错的结果，所以如果是这种情况请告诉我。

ChickWeight[1:20, ]
dim(ChickWeight)
sum(is.na(ChickWeight)) #contains no NAs
ChickWeight$weight[12] <- NA # add 1 NA
ChickWeight$constant <- 1 #add a constant
ChickWeight$Chick <- as.numeric(levels(ChickWeight$Chick)[ChickWeight$Chick]) #class variable has to be an integer

ini <- mice(ChickWeight, maxit = 0)
pred <- ini$predictorMatrix
pred["weight", ] <- c(0, 2, -2, 1, 2)
method <- ini$method
method["weight"] <- "2l.norm"
imputation <- mice(ChickWeight, m = 5, maxit = 5, pred = pred, method = method)

最后一条命令的结果是:

错误 [<-.data.frame ( *tmp* , , i, value = c(37.3233463394145, 159.862324738397 : 替换有 2 行，数据有 1

添加一个额外的 NA 解决了问题

ChickWeight$weight[13] <- NA # add another NA
imputation <- mice(ChickWeight, m = 5, maxit = 5, pred = pred, method = method)

有谁知道可能导致错误的原因是什么？

Answer 1

mice 包的作者 Stef van Buuren 通过电子邮件向我提供了答案。mice 包 v2.22 中的 mice.impute.2l.norm 函数有一个非常小的遗漏，此后已在 v2.23 中修复。 .

正确的代码如下，其中仅在倒数第 4 行添加了 drop = FALSE 到 rowSums(as.matrix(x[nry, type == 2]) , ...:

mice.impute.2l.norm2 <- 
function (y, ry, x, type, intercept = TRUE, ...) 
{
  rwishart <- function(df, p = nrow(SqrtSigma), SqrtSigma = diag(p)) {
    Z <- matrix(0, p, p)
    diag(Z) <- sqrt(rchisq(p, df:(df - p + 1)))
    if (p > 1) {
      pseq <- 1:(p - 1)
      Z[rep(p * pseq, pseq) + unlist(lapply(pseq, seq))] <- rnorm(p * 
                                                                    (p - 1)/2)
    }
    crossprod(Z %*% SqrtSigma)
  }
  force.chol <- function(x, warn = TRUE) {
    z <- 0
    repeat {
      lambda <- 0.1 * z
      XT <- x + diag(x = lambda, nrow = nrow(x))
      XT <- (XT + t(XT))/2
      s <- try(expr = chol(XT), silent = TRUE)
      if (class(s) != "try-error") 
        break
      z <- z + 1
    }
    attr(s, "forced") <- (z > 0)
    if (warn && z > 0) 
      warning("Cholesky decomposition had to be forced", 
              call. = FALSE)
    return(s)
  }
  if (intercept) {
    x <- cbind(1, as.matrix(x))
    type <- c(2, type)
  }
  n.iter <- 100
  nry <- !ry
  n.class <- length(unique(x[, type == (-2)]))
  if (n.class == 0) 
    stop("No class variable")
  gf.full <- factor(x[, type == (-2)], labels = 1:n.class)
  gf <- gf.full[ry]
  XG <- split.data.frame(as.matrix(x[ry, type == 2]), gf)
  X.SS <- lapply(XG, crossprod)
  yg <- split(as.vector(y[ry]), gf)
  n.g <- tabulate(gf)
  n.rc <- ncol(XG[[1]])
  bees <- matrix(0, nrow = n.class, ncol = n.rc)
  ss <- vector(mode = "numeric", length = n.class)
  mu <- rep(0, n.rc)
  inv.psi <- diag(1, n.rc, n.rc)
  inv.sigma2 <- rep(1, n.class)
  sigma2.0 <- 1
  theta <- 1
  for (iter in 1:n.iter) {
    for (class in 1:n.class) {
      vv <- sym(inv.sigma2[class] * X.SS[[class]] + inv.psi)
      bees.var <- chol2inv(chol(vv))
      bees[class, ] <- drop(bees.var %*% (crossprod(inv.sigma2[class] * 
                                                      XG[[class]], yg[[class]]) + inv.psi %*% mu)) + 
        drop(rnorm(n = n.rc) %*% chol(sym(bees.var)))
      ss[class] <- crossprod(yg[[class]] - XG[[class]] %*% 
                               bees[class, ])
    }
    mu <- colMeans(bees) + drop(rnorm(n = n.rc) %*% chol(chol2inv(chol(sym(inv.psi)))/n.class))
    inv.psi <- rwishart(df = n.class - n.rc - 1, SqrtSigma = chol(chol2inv(chol(sym(crossprod(t(t(bees) - 
                                                                                                  mu)))))))
    inv.sigma2 <- rgamma(n.class, n.g/2 + 1/(2 * theta), 
                         scale = 2 * theta/(ss * theta + sigma2.0))
    H <- 1/mean(inv.sigma2)
    sigma2.0 <- rgamma(1, n.class/(2 * theta) + 1, scale = 2 * 
                         theta * H/n.class)
    G <- exp(mean(log(1/inv.sigma2)))
    theta <- 1/rgamma(1, n.class/2 - 1, scale = 2/(n.class * 
                                                     (sigma2.0/H - log(sigma2.0) + log(G) - 1)))
  }
  imps <- rnorm(n = sum(nry), sd = sqrt(1/inv.sigma2[gf.full[nry]])) + 
    rowSums(as.matrix(x[nry, type == 2, drop = FALSE]) * bees[gf.full[nry], 
                                                              ])
  return(imps)
}

当新函数被添加到mice命名空间时

environment(mice.impute.2l.norm2) <- asNamespace('mice')

它可以通过调用 2l.norm2 在小鼠中正常使用。在前面描述的示例中显示:

method["weight"] <- "2l.norm2"
imputation <- mice(ChickWeight, m = 5, maxit = 5, pred = pred, method = method)

它现在按预期工作了!