如何解决绘制polr模型效应图时的报错问题

问题背景 #

使用R语言polr函数构建比例优势逻辑回归模型后，尝试绘制模型效应图时持续报错。数据包含以下变量：

• country：字符型（共5个国家）
• gender：哑变量（1/2）
• age：连续型
• income：1-10缩放的连续型
• 其他变量：Work less、Work much（0/1）、health、education、marital status等

数据样例 #

变量组1 #

变量组2 #

尝试的方法及报错 #

方法1：模型公式内转换因变量类型 #

for.plot <- polr(factor(as.ordered(lifesatisfaction)) ~ country*(gender + age + income + educ + health + `work less` + `work much`), data = surveywave5, method = "logistic", Hess = TRUE)
plot(Effect(focal.predictors = c("country","`work less`"), mod = for.plot, xlevels = list(age = 15:65)), rug = FALSE)

报错信息：

Error in contrasts<-(*tmp*, value = contr.funs[1 + isOF[nn]]): contrasts can be applied only to factors with 2 or more levels

方法2：模型公式内转换分类变量 #

for.plot <- polr(as.factor(as.ordered(lifesatisfaction)) ~ as.factor(country)*(gender + age + income + education + health + as.factor(`work less`) + `work much`), data = surveywave5, method = "logistic", Hess = TRUE)
plot(Effect(focal.predictors = c("country","`work less`"), mod = for.plot, xlevels = list(age = 15:65)), rug = FALSE)

报错信息：

Error in Effect(focal.predictors = c("country", "work less"), mod = for.plot, :
model formula should not contain calls to
factor(), as.factor(), ordered(), as.ordered(), as.numeric(), or as.integer();
see 'Warnings and Limitations' in ?Effect

方法3：提前转换变量类型 #

提前将相关变量转为因子类型，数据结构匹配示例要求，但绘制效应图仍失败。

解决方案 #

步骤1：预处理数据，消除变量类型问题 #

1. 重命名含空格的变量：含空格的变量名会干扰effects包识别，建议重命名为下划线格式：

   library(dplyr)
   surveywave5 <- surveywave5 %>%
     rename(work_less = `Work less`,
            work_much = `Work much`)
   

2. 转换变量类型：

   # 因变量转为有序因子（polr要求因变量为有序因子）
   surveywave5$lifestatisfied <- ordered(surveywave5$lifestatisfied)
   # 分类自变量转为因子
   surveywave5$country <- factor(surveywave5$country)
   surveywave5$gender <- factor(surveywave5$gender)
   surveywave5$work_less <- factor(surveywave5$work_less)
   surveywave5$work_much <- factor(surveywave5$work_much)
   

3. 检查单水平因子：报错“contrasts can be applied only to factors with 2 or more levels”通常意味着存在只有1个水平的因子变量，检查并处理：

   # 查看所有因子变量的水平数
   sapply(surveywave5[, sapply(surveywave5, is.factor)], nlevels)
   # 若存在水平数为1的变量，要么从模型中删除该变量，要么检查数据补全缺失类别
   

4. 处理缺失值：大量NA值可能导致模型或绘图异常，可先删除缺失值或用合适方法插补：

   # 简单删除含NA的行（根据实际情况选择）
   surveywave5_clean <- na.omit(surveywave5)
   

步骤2：重新构建模型（公式中不使用类型转换函数） #

library(MASS)
# 使用清洗后的数据构建模型
for.plot <- polr(lifestatisfied ~ country*(gender + age + income + education + health + work_less + work_much), 
                 data = surveywave5_clean, method = "logistic", Hess = TRUE)

步骤3：绘制效应图 #

library(effects)
# 构建效应对象，聚焦目标变量
eff <- Effect(focal.predictors = c("country", "work_less"), mod = for.plot, xlevels = list(age = 15:65))
# 绘图
plot(eff, rug = FALSE)

额外排查点 #

• 确保effects包为最新版本，执行update.packages("effects")更新
• 若仍报错，可尝试简化模型，先只保留country和work_less两个变量绘图，逐步添加其他变量排查问题

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