brmultinom模型无法运行dredge、VIF及PseudoR2函数的问题求助

我用brmultinom替代原nnet::multinom模型后，原本在multinom上正常运行的dredge（MuMIn包）、vif、PseudoR2（DescTools包）函数均无法正常执行，报错及代码如下：

模型构建代码 #

library(MuMIn)
options(na.action = na.fail)

dat$Sex <- as.factor(dat$Sex)
dat$season <- as.factor(dat$season)
dat$Data.origin <- as.factor(dat$Data.origin)
str(dat$Length)
str(dat$VBT)

# 构建brmultinom模型，prey_new为响应变量矩阵
model2 <- brmultinom(prey_new ~ 
                      Sex +
                      season +
                      Data.origin +
                      Length +
                      VBT,
                    data = dat,type="AS_mean",maxit=999)

1. dredge函数报错 #

执行代码：

summary(model2)
dredge_model_multi <- dredge(model2, rank = "AICc")

错误信息：

Fixed term is "(Intercept)"
Error in cf[, 2L] : subscript out of bounds

2. vif函数报错 #

执行代码：

options(na.action = na.omit)
vif(model2, na.rm = T)

错误信息：

Error in if (names(coefficients(mod)[1]) == "(Intercept)") { :
argument is of length zero

3. PseudoR2函数报错 #

执行代码：

library("DescTools")
# Calculate the R Square
PseudoR2(model2, which = c("CoxSnell","Nagelkerke","McFadden"))

错误信息：

Error in match.arg(type) : 'arg' should be one of “class”, “probs”

解决建议 #

针对dredge函数不兼容问题 #

MuMIn::dredge原生不支持brmultinom模型类，可通过两种方式处理：

• 手动模型选择：遍历所有变量组合，逐个拟合brmultinom模型并计算AICc，自行排序筛选最优模型。
• 自定义S3方法：为brmultinom编写适配dredge的方法，让函数能识别模型结构：

# 自定义logLik方法，提取模型对数似然值
logLik.brmultinom <- function(object, ...) {
  ll <- object$logLik
  attr(ll, "df") <- length(coef(object))
  class(ll) <- "logLik"
  ll
}

# 自定义getAllTerms方法，提取模型公式项
getAllTerms.brmultinom <- function(object, ...) {
  terms(formula(object))
}

添加方法后重新运行dredge即可。

针对vif函数报错问题 #

vif无法识别brmultinom的系数矩阵结构，可手动基于设计矩阵计算VIF：

# 提取不含截距的设计矩阵
X <- model.matrix(prey_new ~ Sex + season + Data.origin + Length + VBT, data = dat)[, -1]

# 计算每个变量的VIF
vif_vals <- apply(X, 2, function(x) {
  # 拟合当前变量对其他所有变量的回归
  fit <- lm(x ~ . - 1, data = as.data.frame(X[, colnames(X) != colnames(X)[which(colnames(X) == deparse(substitute(x)))]])
  r_sq <- summary(fit)$r.squared
  1 / (1 - r_sq)
})

针对PseudoR2函数报错问题 #

DescTools::PseudoR2对brmultinom支持不完善，可手动计算伪R²或换用其他包：

• 手动计算McFadden伪R²：

# 拟合仅含截距的零模型
null_model <- brmultinom(prey_new ~ 1, data = dat, type="AS_mean", maxit=999)
# 提取对数似然值
ll_model <- model2$logLik
ll_null <- null_model$logLik
# 计算McFadden R²
mcfadden_r2 <- 1 - ll_model / ll_null

CoxSnell和Nagelkerke伪R²可参考多项Logistic回归公式，基于对数似然值推导计算。

• 尝试pscl包：

library(pscl)
pR2(model2)

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