如何绘制R中clm()/polr()序数回归模型的交互效应图？

使用ordinal包clm()函数或MASS包polr()函数拟合序数回归模型后，调用interactions包interact_plot()绘制交互效应时会触发报错：

Error in UseMethod("family") : no applicable method for 'family' applied to an object of class 'clm'

换用polr()拟合模型运行相同代码会出现同类报错，复现代码如下：

m <- clm(anxiety_levels ~ pred_a * pred_b + pred_c, data, link = "logit")
interact_plot(m, pred = pred_a, modx = pred_b)

报错核心原因是interact_plot()依赖family()S3方法提取模型分布与链接函数信息，而clm、polr类模型对象没有内置对应方法，无法被函数直接识别。

按稳定性、自定义灵活度从高到低排序：

方案1：基于emmeans计算边际效应后用ggplot2绘图（最推荐） #

该方法适配所有主流序数回归模型，结果可控，不会出现包版本兼容导致的计算错误，步骤如下：

1. 根据调节变量pred_b的类型选典型取值：连续变量取均值、均值±1标准差或常用分位数；分类变量直接取所有观测水平
2. 调用emmeans()计算核心自变量、调节变量不同取值组合下，焦虑各等级的预测概率与置信区间
3. 整理结果后用ggplot2按需绘图
   示例代码：

library(emmeans)
library(ggplot2)

# 设定调节变量pred_b的取值：连续变量用下方写法，分类变量直接传所有水平即可
pred_b_set <- c(
  "低水平(-1SD)" = mean(data$pred_b) - sd(data$pred_b),
  "平均水平" = mean(data$pred_b),
  "高水平(+1SD)" = mean(data$pred_b) + sd(data$pred_b)
)

# 计算各组合下的预测概率
pred_res <- emmeans(
  m,
  specs = ~ pred_a * pred_b,
  at = list(pred_b = pred_b_set),
  mode = "prob" # 直接返回响应变量各等级的预测概率，改为"latent"可输出潜在变量得分
)
pred_df <- as.data.frame(pred_res)

# 绘制交互图，按焦虑等级分面展示
ggplot(pred_df, aes(x = pred_a, y = prob, color = as.factor(pred_b), group = as.factor(pred_b))) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = asymp.LCL, ymax = asymp.UCL, fill = as.factor(pred_b)), alpha = 0.2, color = NA) +
  facet_wrap(~anxiety_levels, nrow = 1, labeller = label_both) +
  labs(x = "核心自变量pred_a", y = "对应焦虑等级的预测概率", color = "调节变量pred_b", fill = "调节变量pred_b") +
  theme_bw()

如果需要展示累积概率、优势比等其他指标，只需要调整emmeans()的参数即可。

方案2：用effects包直接生成交互图 #

effects包原生支持clm、polr类对象，不需要手动计算预测值，代码更简洁，适合快速出图：

library(effects)
# 提取pred_a与pred_b的交互效应
inter_eff <- Effect(focal.predictors = c("pred_a", "pred_b"), mod = m)
# 直接绘图，指定x轴与分组变量即可
plot(
  inter_eff,
  x.var = "pred_a",
  lines.var = "pred_b",
  rug = FALSE,
  ylab = "各焦虑等级预测概率",
  xlab = "pred_a取值"
)

方案3：手动补全S3方法适配interact_plot()（不推荐） #

如果一定要使用interact_plot()，可以临时为clm/polr类补充family方法，让函数可以识别模型类型，但该方法可能因包版本差异出现预测值计算偏差，优先选前两种方案：

# 为clm类添加family方法
family.clm <- function(object, ...) {
  binomial(link = object$link)
}
# 如果用的是polr模型，替换为下方方法
# family.polr <- function(object, ...) {
#   link <- ifelse(object$method == "logistic", "logit", object$method)
#   binomial(link = link)
# }

# 运行绘图函数，注意指定type="prob"输出分类概率
interact_plot(m, pred = pred_a, modx = pred_b, type = "prob")

内容的提问来源于stack exchange，提问作者Rhys Maredudd Davies