使用emmeans的interaction参数时，如何反转对比方向？

问题描述 #

我在零膨胀负二项（ZINB）模型中使用emmeans，通过自定义正交对比检验对比的对比：研究设计包含4组（study_group: grp1、grp2、grp3、grp4），每组在3个时间点（time: Time1、Time2、Time3）评估。

当前代码生成的对比以grp1/grp2、grp1/grp3……grp3/grp4这类「低组/高组」的比率形式呈现，但我需要反转成grp2/grp1、grp3/grp1……grp4/grp3的「高组/低组」形式。尝试在多个位置添加reverse=TRUE均无效，想找除了重新设置study_group因子水平之外的实现方法。

附当前代码及输出：

library(glmmTMB)
library(emmeans)

set.seed(3456)

# 构建研究设计网格：4组，每组3个站点，每个站点20名参与者，观测3次
site <- rep(1:12, each=60)
pid <- 1000*site+10*(rep(rep(1:20,each=3),12))
study_group <- c(rep("grp1",180), rep("grp2",180), rep("grp3",180), rep("grp4",180))
grp_num <- c(rep(0,180), rep(1,180), rep(2,180), rep(3,180))
time <- c(rep(c("Time1", "Time2", "Time3"),240))
time_num <- c(rep(c(0:2),240))

# 站点水平随机效应（截距）
site_eff_count = rep(rnorm(12, mean = 0, sd = 0.5), each = 60)
site_eff_zeros = rep(rnorm(12, mean = 0, sd = 0.5), each = 60)

# 模拟负二项结果
y_count <- rnbinom(n = 720, mu=exp(3.25 + grp_num*0.15 + time_num*-0.20 + grp_num*time_num*0.15 + site_eff_count), size=0.8)

# 模拟额外零值
log_odds = (-1.75 + grp_num*0.2 + time_num*-0.40 + grp_num*time_num*0.50 + site_eff_zeros)
prob_1 = plogis(log_odds)
prob_0 = 1 - prob_1
y_zeros <- rbinom(n = 720, size = 1, prob = prob_0) 

# 构建ZINB数据集
data_ZINB <- data.frame(site, pid, study_group, time, y_count, y_zeros)
data_ZINB$y_obs <- ifelse(y_zeros==1, y_count, 0)

# 拟合ZINB混合模型
mod_ZINB <- glmmTMB(y_obs ~ 1 
                    + study_group + time + study_group*time
                    + (1|site),
                    family=nbinom2,
                    zi = ~ .,
                    data=data_ZINB)

# 获取条件模型（非零部分）的单元格均值（响应尺度）
count_means <- emmeans(mod_ZINB, 
                       pairwise ~ time | study_group, 
                       component="cond", 
                       type="response", 
                       adjust="none")

# 自定义时间正交对比函数：contr1=Time2-Time1；contr2=Time3-(Time1+Time2)/2
compare_arms.emmc <- function(levels) {
  k <- length(levels)
  contr1 <- c(-1,1,0)
  contr2 <- c(-1,-1,2)
  coef <- as.data.frame(lapply(seq_len(k - 1), function(i) {
    if(i==1) contr1 else contr2
  }))
  names(coef) <- c("T1vT2", "T1T2vT3")
  attr(coef, "adjust") = "none"
  coef
}

# 估计组间的「对比的对比」（检验时间对比是否在组间存在差异）
compare_arms_contrast <- contrast(count_means[[1]], 
                                  interaction = c("compare_arms", "pairwise"), 
                                  by = NULL)
compare_arms_contrast

当前输出（截取部分）：

time_compare_arms study_group_pairwise ratio    SE  df null t.ratio p.value
 T1vT2             grp1 / grp2          1.091 0.368 693    1   0.259  0.7957
 T1T2vT3           grp1 / grp2          0.623 0.371 693    1  -0.794  0.4276
 T1vT2             grp1 / grp3          1.190 0.399 693    1   0.520  0.6034
 T1T2vT3           grp1 / grp3          0.384 0.241 693    1  -1.523  0.1283
 ...
 T1T2vT3           grp3 / grp4          0.676 0.556 693    1  -0.475  0.6346

Tests are performed on the log scale

解决方案 #

方法1：自定义反转的组间成对对比函数 #

通过自定义emmc类型的函数，直接生成反转方向的组间对比，替换默认的pairwise函数：

# 自定义反转的组间成对对比函数
reverse_pairwise.emmc <- function(levels) {
    # 获取默认pairwise对比矩阵
    mat <- pairwise.emmc(levels)
    # 反转对比方向（每行取反，对应log比率的符号反转）
    rev_mat <- -mat
    # 修改对比标签，将"A/B"改为"B/A"
    colnames(rev_mat) <- sapply(colnames(rev_mat), function(x) {
        parts <- strsplit(x, " / ")[[1]]
        paste(parts[2], parts[1], sep = " / ")
    })
    # 保留原有的调整方法属性
    attr(rev_mat, "adjust") <- attr(mat, "adjust")
    rev_mat
}

# 使用自定义的反转对比函数生成「对比的对比」
compare_arms_contrast_rev <- contrast(count_means[[1]], 
                                      interaction = c("compare_arms", "reverse_pairwise"), 
                                      by = NULL)
compare_arms_contrast_rev

方法2：事后修改对比结果 #

如果不想修改对比生成逻辑，可以直接对已有的emmGrid结果进行调整：

# 将结果转换为数据框方便修改
contrast_df <- as.data.frame(compare_arms_contrast)

# 反转比率、标签和t值
contrast_df <- contrast_df %>%
    mutate(
        # 反转比率（1/原比率）
        ratio = 1 / ratio,
        # 反转组对比标签
        study_group_pairwise = sapply(study_group_pairwise, function(x) {
            parts <- strsplit(x, " / ")[[1]]
            paste(parts[2], parts[1], sep = " / ")
        }),
        # t值取反（因为log(1/r) = -log(r)，检验逻辑不变，p值不受影响）
        t.ratio = -t.ratio
    ) %>%
    # 保持原有列顺序
    select(time_compare_arms, study_group_pairwise, ratio, SE, df, null, t.ratio, p.value)

# 可选：将修改后的数据重新转为emmGrid对象（保留emmeans后续操作能力）
compare_arms_contrast_rev <- compare_arms_contrast
compare_arms_contrast_rev@grid$study_group_pairwise <- contrast_df$study_group_pairwise
compare_arms_contrast_rev@estimates[, "ratio"] <- contrast_df$ratio
compare_arms_contrast_rev@estimates[, "t.ratio"] <- contrast_df$t.ratio

# 查看反转后的结果
compare_arms_contrast_rev

说明 #

• 两种方法都会得到grp2/grp1、grp3/grp1这类你需要的比率方向，且检验的统计逻辑完全一致（p值不变）。
• 原代码中reverse=TRUE无效的原因：在interaction类型的对比中，reverse参数无法传递给内部的pairwise对比函数，因此需要自定义对比逻辑或事后调整。

内容的提问来源于stack exchange，提问作者M. Todd