如何用emmeans获取mtcars数据集每个观测的disp拟合值？

问题背景 #

我想用emmeans处理mtcars数据集，目前能通过以下代码得到disp在指定值下的平均预测值：

library(emmeans)
mod <- lm(mpg ~ disp + hp, data = mtcars)
emm <- emmeans(mod, "disp", at=list(disp=c(130,135,140)))

但我想知道，能不能像predict(mod)那样，直接获取数据集中每个样本的拟合值？

另外，我尝试了Ben Bolker提供的方法计算均值：

nd <- expand.grid(hp = unique(mtcars$hp), disp = 130)
mean(predict(mod, newdata = nd))
# 结果：23.20037

但这个结果和emmeans的输出不匹配：

emmeans(mod, "disp", at=list(disp=130))

# 输出：
disp emmean    SE df lower.CL upper.CL
  130   23.1 0.928 29     21.2       25

解答1：用emmeans获取每个样本的拟合值 #

完全可以，有两种简单方法：

方法一：基于原数据生成emmeans对象 #

直接指定at参数为原数据集的变量取值，提取emmean就是拟合值：

# 生成每个样本对应的emmeans对象
emm_fitted <- emmeans(mod, ~ disp + hp, at = as.list(mtcars))
# 提取拟合值
fitted_values <- as.data.frame(emm_fitted)$emmean

# 和predict(mod)结果完全一致
all.equal(fitted_values, predict(mod))
# [1] TRUE

方法二：用ref_grid配合predict #

利用emmeans的参考网格（ref_grid），直接调用predict方法：

# 基于原数据创建参考网格
emm_ref <- ref_grid(mod, data = mtcars)
# 获取每个样本的拟合值
predict(emm_ref)

解答2：两种方法结果差异的原因 #

核心是对hp变量的平均权重不同：

• emmeans计算的是边际均值（marginal mean）：它会先对每个hp水平计算预测值，再按照原数据中每个hp的出现次数加权平均。比如mtcars中hp=110出现3次，hp=93只出现1次，emmeans会给前者更高的权重。
• Ben的方法用unique(mtcars$hp)，相当于给每个不同的hp水平赋予了相同权重，完全忽略了原数据中的频数分布，所以计算出的均值和emmeans的边际均值有差异。

验证一下加权平均的结果：

# 按hp的实际频数加权计算
hp_counts <- table(mtcars$hp)
nd_weighted <- expand.grid(hp = names(hp_counts), disp = 130)
preds_weighted <- predict(mod, newdata = nd_weighted)
weighted.mean(preds_weighted, w = as.numeric(hp_counts))
# 结果：23.10357，和emmeans输出的23.1一致

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