一种解决方案是为mice包编写自己的自定义插补函数。该软件包为此做好了准备，并且设置令人惊讶地无痛。
首先，我们按照建议设置数据：
dat=data.frame(x1=c(21, 50, 31, 15, 36, 82, 14, 14, 19, 18, 16, 36, 583, NA,NA,NA, 50, 52, 26, 24),
x2=c(0, NA, 18,0, 19, 0, NA, 0, 0, 0, 0, 0, 0,NA,NA, NA, 22, NA, 0, 0),
x3=c(0, 0, 0, 0, 0, 54, 0 ,0, 0, 0, 0, 0, 0, NA, NA, NA, NA, 0, 0, 0))
接下来我们加载mice包，看看它默认选择了哪些方法：
library(mice)
# Do a non-imputation
imp_base <- mice(dat, m=0, maxit = 0)
# Find the methods that mice chooses
imp_base$method
# Returns: "pmm" "pmm" "pmm"
# Look at the imputation matrix
imp_base$predictorMatrix
# Returns:
# x1 x2 x3
#x1 0 1 1
#x2 1 0 1
#x3 1 1 0
pmm代表预测均值匹配- 可能是用于插补连续变量的最流行的插补算法。它使用回归模型计算预测值，并选择最接近预测值的 5 个元素（通过欧几里得距离）。这些选择的元素称为供体池，最终值是从该供体池中随机选择的。
从预测矩阵中，我们发现这些方法获得了对限制感兴趣的传递变量。请注意，行是目标变量，列是预测变量。如果 x1 在 x3 列中没有 1，我们将不得不将其添加到矩阵中：imp_base$predictorMatrix["x1","x3"] <- 1
现在到有趣的部分，生成插补方法。我在这里选择了一种相当粗略的方法，如果它们不符合标准，我将丢弃所有值。这可能会导致较长的循环时间，并且保留有效的插补并仅重做剩余的插补可能会更有效，尽管这需要更多的调整。
# Generate our custom methods
mice.impute.pmm_x1 <-
function (y, ry, x, donors = 5, type = 1, ridge = 1e-05, version = "",
...)
{
max_sum <- sum(max(x[,"x2"], na.rm=TRUE),
max(x[,"x3"], na.rm=TRUE))
repeat{
vals <- mice.impute.pmm(y, ry, x, donors = 5, type = 1, ridge = 1e-05,
version = "", ...)
if (all(vals < max_sum)){
break
}
}
return(vals)
}
mice.impute.pmm_x2 <-
function (y, ry, x, donors = 5, type = 1, ridge = 1e-05, version = "",
...)
{
repeat{
vals <- mice.impute.pmm(y, ry, x, donors = 5, type = 1, ridge = 1e-05,
version = "", ...)
if (all(vals == 0 | vals >= 14)){
break
}
}
return(vals)
}
mice.impute.pmm_x3 <-
function (y, ry, x, donors = 5, type = 1, ridge = 1e-05, version = "",
...)
{
repeat{
vals <- mice.impute.pmm(y, ry, x, donors = 5, type = 1, ridge = 1e-05,
version = "", ...)
if (all(vals == 0 | vals >= 16)){
break
}
}
return(vals)
}
一旦我们完成了定义方法，我们就可以简单地更改以前的方法。如果您只想更改单个变量，那么您可以简单地使用imp_base$method["x2"] <- "pmm_x2"，但对于此示例，我们将全部更改（命名不是必需的）：
imp_base$method <- c(x1 = "pmm_x1", x2 = "pmm_x2", x3 = "pmm_x3")
# The predictor matrix is not really necessary for this example
# but I use it just to illustrate in case you would like to
# modify it
imp_ds <-
mice(dat,
method = imp_base$method,
predictorMatrix = imp_base$predictorMatrix)
现在让我们看一下第三个估算数据集：
> complete(imp_ds, action = 3)
x1 x2 x3
1 21 0 0
2 50 19 0
3 31 18 0
4 15 0 0
5 36 19 0
6 82 0 54
7 14 0 0
8 14 0 0
9 19 0 0
10 18 0 0
11 16 0 0
12 36 0 0
13 583 0 0
14 50 22 0
15 52 19 0
16 14 0 0
17 50 22 0
18 52 0 0
19 26 0 0
20 24 0 0
好的，这就是工作。我喜欢这个解决方案，因为您可以搭载主流功能，只需添加您认为有意义的限制。
更新
为了强制执行评论中提到的@t0x1n 的严格限制，我们可能希望在包装函数中添加以下功能：
- 在循环期间保存有效值，以便之前部分成功运行的数据不会被丢弃
- 一种逃避机制，以避免无限循环
- 在尝试x次但没有找到合适的匹配项后膨胀供体池（这主要适用于 pmm）
这导致了一个稍微复杂的包装函数：
mice.impute.pmm_x1_adv <- function (y, ry,
x, donors = 5,
type = 1, ridge = 1e-05,
version = "", ...) {
# The mice:::remove.lindep may remove the parts required for
# the test - in those cases we should escape the test
if (!all(c("x2", "x3") %in% colnames(x))){
warning("Could not enforce pmm_x1 due to missing column(s):",
c("x2", "x3")[!c("x2", "x3") %in% colnames(x)])
return(mice.impute.pmm(y, ry, x, donors = 5, type = 1, ridge = 1e-05,
version = "", ...))
}
# Select those missing
max_vals <- rowSums(x[!ry, c("x2", "x3")])
# We will keep saving the valid values in the valid_vals
valid_vals <- rep(NA, length.out = sum(!ry))
# We need a counter in order to avoid an eternal loop
# and for inflating the donor pool if no match is found
cntr <- 0
repeat{
# We should be prepared to increase the donor pool, otherwise
# the criteria may become imposs
donor_inflation <- floor(cntr/10)
vals <- mice.impute.pmm(y, ry, x,
donors = min(5 + donor_inflation, sum(ry)),
type = 1, ridge = 1e-05,
version = "", ...)
# Our criteria check
correct <- vals < max_vals
if (all(!is.na(valid_vals) |
correct)){
valid_vals[correct] <-
vals[correct]
break
}else if (any(is.na(valid_vals) &
correct)){
# Save the new valid values
valid_vals[correct] <-
vals[correct]
}
# An emergency exit to avoid endless loop
cntr <- cntr + 1
if (cntr > 200){
warning("Could not completely enforce constraints for ",
sum(is.na(valid_vals)),
" out of ",
length(valid_vals),
" missing elements")
if (all(is.na(valid_vals))){
valid_vals <- vals
}else{
valid_vals[is.na(valid_vals)] <-
vals[is.na(valid_vals)]
}
break
}
}
return(valid_vals)
}
请注意，这并没有表现得那么好，很可能是因为建议的数据集在所有情况下都没有通过约束而没有丢失。我需要在它开始表现之前将循环长度增加到 400-500。我认为这是无意的，您的估算应该模仿实际数据的生成方式。
优化
该参数ry包含非缺失值，我们可以通过删除我们发现符合条件的插补的元素来加速循环，但由于我不熟悉内部函数，所以我避免这样做。
我认为当您有需要时间来完成的强大约束时，最重要的事情是并行化您的插补（请参阅我在 CrossValidated 上的回答）。大多数今天的计算机都有 4-8 个内核，而 R 默认只使用其中一个。通过将内核数量增加一倍，时间可以（几乎）减半。
插补时缺少参数
关于在x2插补时丢失的问题 - 老鼠实际上从不将缺失值输入x- data.frame。老鼠方法包括在开始时填写一些随机值。插补的链部分限制了该初始值的影响。如果您查看-function，您可以在插补调用（ -function）mice之前找到它：mice:::sampler
...
if (method[j] != "") {
for (i in 1:m) {
if (nmis[j] < nrow(data)) {
if (is.null(data.init)) {
imp[[j]][, i] <- mice.impute.sample(y,
ry, ...)
}
else {
imp[[j]][, i] <- data.init[!ry, j]
}
}
else imp[[j]][, i] <- rnorm(nrow(data))
}
}
...
data.init可以提供给mice函数，mice.imput.sample是一个基本的采样过程。
参观顺序
如果访问顺序很重要，您可以指定mice-function 运行插补的顺序。默认为 from1:ncol(data)但您可以将 设置为visitSequence您喜欢的任何内容。