如何解决rmultinom报错：概率向量中存在NA？

问题重现 #

使用VIM包的colic数据集复现报错，代码如下：

library(VIM)

data(colic)
summary(colic)

###########subset data
newcolic <- colic[,-c(1,2,3,7,12,13,14,15,16,17,19,20,21,24,25,26,27, 29), drop=FALSE] 
summary(newcolic)
str(newcolic)

############### reorder categories
levels(newcolic$pulse_peripheral)
newcolic$pulse_peripheral <- factor(newcolic$pulse_peripheral, levels = c("normal", "increased", "reduced", "absent"))
summary(newcolic$pulse_peripheral)

newcolic$peristalsis <- as.factor(newcolic$peristalsis)
is.factor(newcolic$peristalsis)

str(newcolic)

#########  Combine variables
library(forcats)
library(dplyr)
newcolic <- (newcolic %>% mutate(outcome = fct_collapse(outcome, "dead" = c("died", "euthanized"))))
summary(newcolic$outcome)
levels(newcolic$outcome)

################################# rename levels
levels(newcolic$mucous_membranes_group)[levels(newcolic$mucous_membranes_group)=="serios"] <- "serious"
levels(newcolic$mucous_membranes_group)
summary(newcolic$mucous_membranes_group)

summary(newcolic)

####Test smcfcs
library(smcfcs)
set.seed(123)

#need to specify method for every variable in dataset
meth <- c("norm", "norm", "norm", "mlogit", "mlogit", "norm", "mlogit", "norm", "", "", "mlogit", "mlogit", "mlogit")

#Error if outcome variable is not coded as numeric
summary(newcolic$outcome)
is.numeric(newcolic$outcome)
newcolic$outcome <- as.numeric(newcolic$outcome)
newcolic$outcome[newcolic$outcome == 2] <- 0

class(newcolic$outcome)
summary(newcolic$outcome)

library(mice)
predictors <- quickpred(newcolic)

#model
imps <- smcfcs(newcolic, smtype="logistic", smformula="outcome ~ temp_rectal + pulse_peripheral + capillayr_refill_time + pain + peristalsis + surgical_lesion + temp_extreme_ordered + mucous_membranes_group",
               method=meth, predictorMatrix=predictors, m=20, numit=15, rjlimit=1000)

运行后报错：

Error in rmultinom(1, size = 1, prob = prob) : NA in probability vector

Traceback信息：

5: rmultinom(1, size = 1, prob = prob)
4: FUN(newX[, i], ...)
3: apply(directImpProbs, 1, catdraw)
2: smcfcs.core(originaldata, smtype, smformula, method, predictorMatrix, 
       m, numit, rjlimit, noisy, errorProneMatrix = errorProneMatrix, 
       restrictions = restrictions)
1: smcfcs(newcolic, smtype = "logistic", smformula = "outcome ~ temp_rectal + pulse_peripheral + capillayr_refill_time + pain + peristalsis + surgical_lesion + temp_extreme_ordered + mucous_membranes_group", 
       method = meth, predictorMatrix = predictors, m = 20, numit = 15, 
       rjlimit = 1000)

报错原因 #

1. 变量拼写错误：模型公式中的capillayr_refill_time是拼写错误，正确变量名为capillary_refill_time，导致模型无法识别该变量，拟合时生成NA概率值。
2. 自动生成的预测矩阵存在无效关系：quickpred生成的predictorMatrix可能包含不合理的预测变量组合，比如让结局变量outcome作为其他变量的预测因子，或在多分类插补模型中引入导致完全分离的变量，使得回归模型无法输出有效概率。
3. 插补方法与变量类型不匹配：若meth参数的元素数量、对应变量的插补方法与newcolic的变量不匹配，会导致部分变量插补逻辑错误，生成NA概率。
4. 分类变量稀有水平：部分分类变量的水平样本量极小，在插补迭代过程中出现完全预测的情况，导致多分类回归拟合失败，输出NA概率。

解决方案 #

1. 修正变量拼写错误 #

将模型公式中的capillayr_refill_time改为正确的capillary_refill_time：

smformula="outcome ~ temp_rectal + pulse_peripheral + capillary_refill_time + pain + peristalsis + surgical_lesion + temp_extreme_ordered + mucous_membranes_group"

2. 手动构建预测矩阵 #

避免使用quickpred自动生成，手动指定合理的预测关系，比如禁止结局变量outcome作为其他变量的预测因子：

# 初始化预测矩阵
predictors <- matrix(1, ncol=ncol(newcolic), nrow=ncol(newcolic))
rownames(predictors) <- colnames(newcolic)
colnames(predictors) <- colnames(newcolic)

# 结局变量不作为其他变量的预测因子
predictors[, "outcome"] <- 0
# 变量不预测自身
diag(predictors) <- 0

3. 核对插补方法与变量的匹配性 #

逐个核对newcolic的变量类型，确保meth参数对应正确：

• 连续变量使用"norm"
• 多分类变量使用"mlogit"
• 结局变量outcome对应的meth元素设为空字符串""

4. 处理分类变量稀有水平 #

检查分类变量的水平分布，合并或删除样本量极小的水平，比如：

# 示例：合并pulse_peripheral的稀有水平
newcolic$pulse_peripheral <- fct_lump_min(newcolic$pulse_peripheral, min=5)

5. 调整拒绝采样限制 #

适当增大rjlimit参数，减少拒绝采样过程中的异常：

imps <- smcfcs(newcolic, 
               smtype="logistic", 
               smformula="outcome ~ temp_rectal + pulse_peripheral + capillary_refill_time + pain + peristalsis + surgical_lesion + temp_extreme_ordered + mucous_membranes_group",
               method=meth, 
               predictorMatrix=predictors, 
               m=20, 
               numit=15, 
               rjlimit=2000)

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