R中polr建模提示设计秩不足丢弃系数如何保留所有变量

运行MASS包的polr函数拟合有序逻辑回归时，出现秩缺陷警告，代码与警告信息如下：

mROBOT_APPRAISAL <- polr(robot_appraisal ~ genderf + up_to_15 + still_studying
                + self_employed + manager + other_whiteC 
                + manual + house_person + unempl + retired,
                 Hess = TRUE,  na.action = na.exclude )
# 警告信息
# Warning message:
# In polr(robot_appraisal ~ genderf + up_to_15 + still_studying +  :
#   design appears to be rank-deficient, so dropping some coefs

默认情况下polr会自动丢弃共线的变量，summary()输出时不会展示被丢弃变量的结果，需要调整处理才能展示所有纳入的变量。

设计矩阵秩缺陷本质是完全多重共线性：即某个变量可以被其他变量的线性组合完全预测，这类变量不存在独立的系数估计值。你当前模型的共线问题基本是两种情况：

• 手动生成的互斥分类哑变量全部放入模型，和截距项共同构成了“虚拟变量陷阱”：比如就业状态是互斥分类，你把所有分类的0-1哑变量都放进模型，所有哑变量的和对于每个样本都是1，和截距项完全共线。
• 某个分类变量在样本中没有变异：比如样本里没有退休状态的受访者，retired变量全为0，无法估计系数。

方案1：保留原模型设定，让summary输出所有变量（共线项显示NA） #

如果不需要调整模型结构，只是希望输出时列出所有你写进公式的变量，可以手动给被丢弃的共线变量补充NA值，修改模型对象的系数存储位置即可：

# 先运行原模型
mROBOT_APPRAISAL <- polr(robot_appraisal ~ genderf + up_to_15 + still_studying
                + self_employed + manager + other_whiteC 
                + manual + house_person + unempl + retired,
                 Hess = TRUE,  na.action = na.exclude )

# 先定位被删掉的共线变量
X <- model.matrix(mROBOT_APPRAISAL)
dropped_vars <- colnames(X)[qr(X)$pivot[(qr(X)$rank+1):ncol(X)]]

# 生成包含所有变量的系数向量，共线项赋值为NA
full_coefs <- setNames(rep(NA, length(c(colnames(X)[-1]))), colnames(X)[-1])
full_coefs[names(coef(mROBOT_APPRAISAL))] <- coef(mROBOT_APPRAISAL)

# 替换模型里的系数对象
mROBOT_APPRAISAL$coefficients <- full_coefs

# 此时调用summary就会展示所有变量，共线项会显示NA，标记为aliased
summary(mROBOT_APPRAISAL)

方案2：从根源解决共线问题，输出有效估计结果 #

如果需要得到所有可估计变量的有效结果，建议不要手动生成互斥哑变量，直接把多分类变量转换为因子类型放入模型，R会自动选择一个分类作为参照组，既避免虚拟变量陷阱，输出也会清晰标注参照项：

# 把互斥的就业/职业哑变量合并为一个因子变量
# 先给occ变量赋值，每个样本匹配对应的职业状态
your_data$occ <- NA
your_data$occ[your_data$still_studying == 1] <- "still_studying"
your_data$occ[your_data$self_employed == 1] <- "self_employed"
your_data$occ[your_data$manager == 1] <- "manager"
your_data$occ[your_data$other_whiteC == 1] <- "other_whiteC"
your_data$occ[your_data$manual == 1] <- "manual"
your_data$occ[your_data$house_person == 1] <- "house_person"
your_data$occ[your_data$unempl == 1] <- "unempl"
your_data$occ[your_data$retired == 1] <- "retired"
# 转为因子，可通过调整levels顺序更换参照组
your_data$occ <- factor(your_data$occ, levels = c("你要设置为参照组的分类名",
                                                  "still_studying","self_employed","manager",
                                                  "other_whiteC","manual","house_person",
                                                  "unempl","retired"))

# 重新拟合模型，不会再出现秩缺陷警告
mROBOT_APPRAISAL_fixed <- polr(robot_appraisal ~ genderf + up_to_15 + occ,
                               Hess = TRUE, na.action = na.exclude, data = your_data)
summary(mROBOT_APPRAISAL_fixed)

注意：完全共线的变量不存在可识别的系数估计值，哪怕强制在输出中展示，对应项也只能显示为NA，不要强行给共线项赋值估计，会导致模型结果完全失效。如果是零单元格导致的共线，需要检查样本分布，考虑合并样本量过少的分类。

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