如何在emmeans的contrast函数中对逆变换结果进行回转换？

问题场景 #

使用faraway包的rats数据集，拟合响应变量经逆变换的线性模型：

data("rats", package = "faraway")
m <- lm(1/time ~ treat + poison, data = rats)

拟合后，能成功获取经回转换的响应尺度边际均值：

emmeans(m, specs = "treat", type = "response", tran = "inverse")

输出结果：

treat response     SE df lower.CL upper.CL
 A        0.284 0.0115 42    0.263    0.309
 B        0.537 0.0411 42    0.465    0.635
 C        0.339 0.0164 42    0.309    0.376
 D        0.463 0.0305 42    0.408    0.534

Results are averaged over the levels of: poison 
Confidence level used: 0.95 
Intervals are back-transformed from the inverse scale 

但计算成对对比时，结果始终停留在逆变换尺度（1/time尺度）：

emmeans(m, specs = "treat", type = "response", tran = "inverse") |> 
  contrast("pairwise")

输出结果：

contrast estimate    SE df t.ratio p.value
 A - B       1.657 0.201 42   8.233  <.0001
 A - C       0.572 0.201 42   2.842  0.0335
 A - D       1.358 0.201 42   6.747  <.0001
 B - C      -1.085 0.201 42  -5.391  <.0001
 B - D      -0.299 0.201 42  -1.485  0.4551
 C - D       0.786 0.201 42   3.905  0.0018

Results are averaged over the levels of: poison 
Note: contrasts are still on the inverse scale 
P value adjustment: tukey method for comparing a family of 4 estimates 

尝试将type = "response", tran = "inverse"移至contrast()函数，得到异常结果：

emmeans(m, specs = "treat") |> 
  contrast("pairwise", type = "response", tran = "inverse")

输出结果：

contrast response  SE df null t.ratio p.value
 A - B           1   0 42  Inf   8.233  <.0001
 A - C           2   1 42  Inf   2.842  0.0335
 A - D           1   0 42  Inf   6.747  <.0001
 B - C         Inf Inf 42  Inf  -5.391  <.0001
 B - D         Inf Inf 42  Inf  -1.485  0.4551
 C - D           1   0 42  Inf   3.905  0.0018

使用make.tran定义变换后拟合模型，对比结果仍处于逆变换尺度：

tran <- make.tran("inverse")
m2 <- with(tran, lm(linkfun(time) ~ treat + poison, data = rats))
emmeans(m2, specs = "treat", type = "response") |> 
  contrast("pairwise")

解决方案 #

要获取原始响应尺度（time尺度）的成对对比结果，需明确：直接逆变换1/time尺度的对比结果是不合理的，因为E[1/time] ≠ 1/E[time]。正确的做法是基于响应尺度的边际均值计算对比，可通过以下两种方式实现：

方法1：直接拟合响应尺度的线性模型（若假设成立） #

如果time的残差满足线性模型的正态性、方差齐性假设，直接在原始尺度拟合模型后计算对比：

m_resp <- lm(time ~ treat + poison, data = rats)
emmeans(m_resp, specs = "treat") |> contrast("pairwise", type = "response")

方法2：转换边际均值后计算对比 #

若必须使用逆变换模型，可先获取响应尺度的边际均值，再用regrid()函数转换为响应尺度的网格对象，最后计算对比：

# 获取响应尺度的边际均值
emm_resp <- emmeans(m, specs = "treat", type = "response", tran = "inverse")
# 转换为响应尺度的网格，再做成对对比
emm_regrid <- regrid(emm_resp)
contrast(emm_regrid, "pairwise")

该方法会基于原始time尺度的边际均值计算成对差异，同时正确推导标准误和置信区间，结果符合预期。

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