finalfit包or_plot函数无法同时纳入主效应与交互项的问题求助

问题现象 #

使用finalfit的or_plot函数时，遇到以下异常：

• 仅包含主效应的模型：可正常生成森林图
• 仅包含交互项（不含对应主效应）的模型：可正常运行
• 同时包含主效应与交互项的模型：要么交互项不显示，要么触发*dplyr::tibble()列尺寸不兼容*错误

测试代码 #

library(finalfit)

set.seed(2025)
n <- 100

# 生成分类预测变量
sex <- factor(sample(c("Male", "Female"), n, replace = TRUE))
smoking <- factor(sample(c("Never", "Former", "Current"), n, replace = TRUE))
region <- factor(sample(c("North", "South", "East", "West"), n, replace = TRUE))

# 模拟logistic回归结局
logit_p <- -2 + 
  ifelse(sex == "Male", 0.8, 0) +
  ifelse(smoking == "Current", 1.2, ifelse(smoking == "Former", 0.6, 0)) +
  ifelse(region == "South", 0.5, ifelse(region == "West", 0.3, 0))

p <- 1 / (1 + exp(-logit_p))
outcome <- rbinom(n, 1, prob = p)

# 合并数据集
df <- data.frame(sex, smoking, region, outcome)

# 仅主效应模型：正常运行
explanatory <- c("sex", "smoking")
dependent <- "outcome"
multivarplot <- df %>% or_plot(dependent, explanatory)
plot(multivarplot)

# 仅交互项模型：正常运行
df$sex_smoking_interaction <- paste(df$sex, df$smoking)
explanatory <- c("sex_smoking_interaction", "region")
model_works <- glm(outcome ~ sex_smoking_interaction + region, data = df, family = binomial )
multivarplot_works <- df %>% or_plot(dependent, explanatory)
plot(multivarplot_works)

# 主效应+交互项：交互项不显示
explanatory <- c("sex:smoking", "sex", "smoking")
model_no_display <- glm(outcome ~ sex:smoking + smoking + sex, data = df, family = binomial )
multivarplot_no_display <- df %>% or_plot(dependent, explanatory)

# 主效应+交互项：触发列尺寸错误
explanatory <- c("sex_smoking_interaction", "sex", "smoking")
model_error <- glm(outcome ~ sex_smoking_interaction + smoking + sex, data = df, family = binomial )
multivarplot_error <- df %>% or_plot(dependent, explanatory)

解决方案 #

1. finalfit框架内的修复方法 #

or_plot通过explanatory参数指定变量时，无法正确处理主效应与交互项的维度匹配问题。改用直接传入拟合好的模型对象即可解决：

# 拟合包含主效应+交互项的完整模型（用*自动包含主效应和交互项）
model_full <- glm(outcome ~ sex * smoking + region, data = df, family = binomial)

# 直接将模型对象传入or_plot，强制包含交互项
plot_full <- or_plot(model_full, include_interaction = TRUE)
plot(plot_full)

2. 替代方案：使用forestplot包 #

如果需要更灵活的交互项展示，推荐使用forestplot包，它对复杂模型结果的支持更友好：

library(forestplot)
library(broom)

# 提取模型的OR值、置信区间
tidy_result <- tidy(model_full, exponentiate = TRUE, conf.int = TRUE)

# 自定义变量标签（根据模型系数名称调整）
var_labels <- c(
  "性别（Male vs Female）",
  "吸烟状态（Former vs Never）",
  "吸烟状态（Current vs Never）",
  "地区（South vs North）",
  "地区（East vs North）",
  "地区（West vs North）",
  "性别×吸烟（Male:Former vs 参照组）",
  "性别×吸烟（Male:Current vs 参照组）"
)

# 绘制森林图
forestplot(
  labeltext = var_labels,
  mean = tidy_result$estimate,
  lower = tidy_result$conf.low,
  upper = tidy_result$conf.high,
  xlog = TRUE, # OR值使用对数坐标轴
  title = "Logistic回归结果（含主效应与交互项）",
  xticks = c(0.25, 0.5, 1, 2, 4)
)

3. 替代方案：ggplot2手动绘制 #

如果需要高度自定义的可视化效果，用ggplot2直接构建森林图：

library(ggplot2)
library(broom)

# 提取模型结果并整理
tidy_result <- tidy(model_full, exponentiate = TRUE, conf.int = TRUE)
tidy_result$term <- factor(tidy_result$term, levels = rev(tidy_result$term)) # 反转变量顺序

# 绘制森林图
ggplot(tidy_result, aes(x = estimate, y = term)) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "gray50") +
  geom_point(size = 3, color = "#2c3e50") +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2, color = "#3498db") +
  scale_x_log10(breaks = c(0.25, 0.5, 1, 2, 4)) +
  labs(x = "优势比（OR）", y = "变量", title = "含主效应与交互项的回归森林图") +
  theme_minimal() +
  theme(panel.grid.major.y = element_blank())

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