使用interact_plot绘制三向交互图报错：仅多水平因子可应用对比

问题场景 #

使用interactions包的interact_plot函数基于svyglm对象绘制三向交互图时，未添加额外因子协变量时可正常运行，添加后报错contrasts can be applied only to factors with 2 or more levels。

可复现代码 #

先加载依赖与数据：

library(survey)
library(interactions)
data(api) # 加载内置数据集api

构建调查设计与基础模型（可正常生成交互图）：

# 构建调查对象
dclus2 <- svydesign(id=~dnum+snum, weights=~pw, data=apiclus2)
rclus2 <- as.svrepdesign(dclus2)

# 带三向交互的svyglm模型
l <- svyglm(snum ~ stype * sch.wide * api99, family= gaussian, rclus2)

# 生成交互图（正常运行）
plot1 <- interact_plot(l, pred = api99 , modx = stype, mod2 = sch.wide)

添加因子协变量comp.imp后报错：

l2 <- svyglm(snum ~ stype * sch.wide * api99 + comp.imp, family= gaussian, rclus2)
interact_plot(l2, pred = api99 , modx = stype, mod2 = sch.wide)
# 报错信息：
# Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
#   contrasts can be applied only to factors with 2 or more levels

替换为数值协变量api00则正常运行：

l3 <- svyglm(snum ~ stype * sch.wide * api99 + api00, family= gaussian, rclus2)
interact_plot(l3, pred = api99 , modx = stype, mod2 = sch.wide) # 正常生成图

解决方案 #

问题核心是interact_plot处理svyglm对象时，会自动将协变量固定在典型值（如因子参考水平、数值变量均值），但调查数据的抽样权重可能导致因子在计算典型值的子集内只剩单一水平，触发对比矩阵错误。

方法1：手动指定协变量固定值 #

用covariate.values参数明确指定因子协变量的取值，避免自动计算时出现单水平问题：

# 查看comp.imp的水平
levels(apiclus2$comp.imp) # 输出: "No" "Yes"

# 手动指定comp.imp固定为某一水平（比如"No"）
interact_plot(l2, pred = api99, modx = stype, mod2 = sch.wide,
              covariate.values = list(comp.imp = "No"))

方法2：预处理因子协变量确保水平完整 #

强制因子保留全数据集的所有水平，避免加权后丢失水平：

# 强制comp.imp为因子并保留所有水平
apiclus2$comp.imp <- factor(apiclus2$comp.imp, levels = c("No", "Yes"))

# 重新构建调查设计和模型
dclus2 <- svydesign(id=~dnum+snum, weights=~pw, data=apiclus2)
rclus2 <- as.svrepdesign(dclus2)
l2 <- svyglm(snum ~ stype * sch.wide * api99 + comp.imp, family= gaussian, rclus2)

# 再调用interact_plot
interact_plot(l2, pred = api99 , modx = stype, mod2 = sch.wide)

方法3：用emmeans计算边际效应后绘图 #

绕开interact_plot的自动处理逻辑，手动计算边际效应再可视化：

library(emmeans)

# 计算边际效应，指定要展示的变量组合，固定协变量值
emm <- emmeans(l2, ~ api99 * stype * sch.wide, 
               at = list(api99 = quantile(apiclus2$api99, c(0.1, 0.5, 0.9))),
               cov.reduce = list(comp.imp = "No"))

# 转换为数据框后用ggplot绘图
emm_df <- as.data.frame(emm)
library(ggplot2)
ggplot(emm_df, aes(x = api99, y = emmean, color = stype, linetype = sch.wide)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = stype, color = NA), alpha = 0.2) +
  theme_minimal()

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