如何在ANCOVA的emmeans检验中纳入分类协变量Gender与连续协变量Age？

一、先明确变量类型与模型逻辑 #

你的模型Feature ~ Age + Gender + Category中，Age是连续协变量，而Gender虽然以0/1哑变量形式存储，但本质是二分类分类变量，并非连续协变量。如果直接以数值型0/1输入，R会误将其当作连续变量处理，这是导致emmeans_test报错的核心原因之一。

二、正确纳入双协变量的emmeans_test用法 #

1. 修正Gender的变量类型 #

先将Gender转换为因子，明确其分类属性：

# 转换为因子，指定水平和标签（可选，便于结果解释）
df_copy$Gender <- factor(df_copy$Gender, levels = c(0, 1), labels = c("Male", "Female"))

2. 重新拟合ANCOVA模型 #

res.aov <- df_copy %>% anova_test(Feature ~ Age + Gender + Category)

3. 用emmeans_test同时控制Age和Gender #

如果目标是检验不同Category下Feature的差异，同时控制Age均值与Gender的不同水平，有两种实用方式：

方式1：按Gender分组统一控制Age均值

# 检验Category的两两差异，控制Age均值，同时按Gender分组计算
emm_res <- df_copy %>%
  emmeans_test(
    Feature ~ Category,
    covar = Age,  # 连续协变量指定在covar
    model = res.aov,
    p.adjust.method = "bonferroni",
    by = Gender  # 按Gender拆分结果
  )

方式2：指定Gender的特定水平单独检验

如果需要分别在Male或Female群体中分析Category的差异，用at参数锁定Gender水平：

# 在Male群体中，控制Age均值，检验Category差异
emm_res_male <- df_copy %>%
  emmeans_test(
    Feature ~ Category,
    covar = Age,
    at = list(Gender = "Male"),
    model = res.aov,
    p.adjust.method = "bonferroni"
  )

# 在Female群体中，控制Age均值，检验Category差异
emm_res_female <- df_copy %>%
  emmeans_test(
    Feature ~ Category,
    covar = Age,
    at = list(Gender = "Female"),
    model = res.aov,
    p.adjust.method = "bonferroni"
  )

三、哑变量Gender的额外处理要点 #

1. 必须转换为因子类型：0/1数值型哑变量会被R当作连续协变量，导致边际均值计算在无意义的中间值（如0.5），转换为因子后才能按分类水平正确计算。
2. 明确参考水平：转换因子时指定levels参数，确保参考组符合你的研究逻辑（比如以Male为基准）。
3. 避免直接放入covar参数：Gender作为分类变量，单独放入covar会计算所有Gender的平均水平（无实际意义），建议用by分组或at指定具体水平。

四、备选方案：直接使用emmeans包（更灵活） #

如果ggpubr的emmeans_test封装不够灵活，可直接调用底层的emmeans包：

library(emmeans)

# 计算不同Category在各Gender水平、Age均值下的边际均值
emm <- emmeans(res.aov, ~ Category | Gender, cov.reduce = mean)

# 进行两两比较并调整p值
pairs(emm, adjust = "bonferroni")

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