在R的lm()中实现分组内因子水平均值对比的方法问询

我有一组数据，想要对比某因子的一个水平与该因子另一水平的均值，同时需要在另一因子的所有水平下完成这个对比——也就是获取第二个因子各水平内，关键因子两个水平的均值差异系数。理想状态下希望通过R中单次调用lm()实现该需求，不想依赖边际效应包或单独跑回归。

我尝试过x + x:y的模型语法，单个交互项时能得到预期结果，但加入另一组与cyl的交互项后，模型估计值和边际对比结果出现了偏离。我清楚该现象的原因，但不确定能否通过因子重编码或对比编码解决这个问题。

我的目标是估计am各水平下，cyl=8与cyl=4的均值差异，示例代码如下：

library(marginaleffects) # 用于计算边际对比
data(mtcars)

mtcars$cyl <- factor(mtcars$cyl)
mtcars$gear <- factor(mtcars$gear)
mtcars$am <- factor(mtcars$am)

# 这个模型的回归结果显然不是我要的边际效应
mod <- lm(mpg ~ cyl * am + gear, data = mtcars)
summary(mod)

# 这是我想直接在模型中估计的目标结果：
avg_comparisons(mod,
  variables = "cyl",
  by = "am")

# 这个模型的回归系数刚好能给出我想要的估计值：
mod2 <- lm(mpg ~ am / cyl + gear, data = mtcars)
summary(mod2)

# 对比验证
avg_comparisons(mod2,
  variables = "cyl",
  by = "am")

# 但加入另一组与cyl的交互项后，估计值和边际对比结果就偏离了：
mod3 <- lm(mpg ~  am / cyl + gear / cyl, data = mtcars)
summary(mod3)

avg_comparisons(mod3,
  variables = "cyl",
  by = "am")

# 我找不到修正这个模型的方法

核心问题 #

当模型中同时包含多组与关键因子（cyl）的交互项时，如何通过调整因子编码或模型语法，让lm()的回归系数直接对应目标边际差异（am各水平下cyl=8与cyl=4的均值差）？

内容的提问来源于stack exchange，提问作者T-Porsch