为何特定自变量序列会导致R的GLM模型出现系数输出bug？

核心原因 #

你遇到的现象不是glm.nb函数的异常，是变量构造时的操作失误导致的：

• 构造d向量时，你将缺失值写为了字符串"NA"，而非R可识别的原生缺失值NA，因此e转为因子时会被识别为3水平因子，三个水平分别为"NA"、"No"、"Yes"，并非你预期的二分类因子。
• R的因子默认按字符串字母顺序设置参考水平，"NA"的字母排序优先级高于"No"，因此模型默认将"NA"作为参考水平，自然会同时输出剩余两个水平"No"和"Yes"的系数，和你预期的以"No"为参考的二分类编码结果不符。

验证方法 #

在转因子后运行以下代码即可确认变量实际水平：

# 查看所有因子水平
levels(df$e)
# 查看因子水平总数
nlevels(df$e)

修复方案 #

方案1：构造数据时使用原生缺失值 #

将d向量的构造代码修改为如下形式，使用不加引号的NA表示缺失值：

d <- c("Yes","No","Yes","No",NA,"Yes")

此时转因子后R会自动将真正的NA识别为缺失值，不会纳入因子水平，e即为二分类因子，默认参考水平为"No"，仅会输出"Yes"的系数。

方案2：对已有数据做清洗 #

如果已经生成了带"NA"字符串的数据集，可以通过以下代码修正：

# 把字符串"NA"转为R可识别的原生缺失值
df$e[df$e == "NA"] <- NA
# 清除因子中无观测的多余水平
df$e <- droplevels(df$e)
# 手动指定参考水平为"No"，符合你的分析预期
df$e <- relevel(df$e, ref = "No")

修复后重新拟合模型即可得到符合预期的哑变量编码结果。

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