无法通过 R 数据代码输出

编程问答 2022-05-04

问题描述

我想测试观察者对不同极性（正（白色）和负（黑色）对比刺激的反应）的对称性。我将反应时间 (RT) 作为因变量，有四种不同的对比。众所周知，响应时间遵循 Pieron 曲线，其渐近线位于 (1) 观察者阈值 (Inf) 和 (2) 位于 250 和 450 毫秒之间的基本 RT 处。 该知识使我们能够将自变量（有效对比度 EC）的关系线性化为 1/EC^2 (tEC)，因此将 RT 与 EC 联系起来的方程变为：

RT = m * tEC + RT0

为了测试对称性，我建立了标准：两个极性中的相同斜率和相同截距意味着对称。 为了获得系数，我制作了一个带有交互作用的线性模型（通过一个虚拟变量编码极性：正或负）。 lm 的输出对我来说很清楚，但有些同事更喜欢与方差分析输出更相似的东西。所以我决定用emmeans来做对比。斜率没问题，但是当计算截距时，问题就开始了。 lm计算的截距和emmeans的输出有很大的不同，结论也不同。在接下来的内容中，我重现了这个例子。 问题有两个：是否可以使用emmeans来解决我的问题？如果没有，是否可以通过其他包（哪个）进行对比？

数据

RT1000 EC tEC 极性
596.3564 -25 0.001600 否定
648.2471 -20 0.002500 否定
770.7602 -17 0.003460 否定
831.2971 -15 0.004444 否定
1311.3331 15 0.004444 正面
1173.8942 17 0.003460 正面
1113.7240 20 0.002500 正面
869.3635 25 0.001600 正面

代码

# Model
model <- lm(RT1000 ~ tEC * Polarity,data = Data)

# emmeans
library(emmeans)
# Slopes
m.slopes <- lstrends(model,"Polarity",var="tEC")
# Intercepts
m.intercept <- lsmeans(model,"Polarity")

# Contrasts
pairs(m.slopes)
pairs(m.intercept)

输出

模型
任期 估计 标准错误 统计 p.value
（拦截） 449.948 66.829 6.733 0.003
tEC 87205.179 20992.976 4.154 0.014
PolarityPositive 230.946 94.511 2.444 0.071
tEC:PolarityPositive 58133.172 29688.551 1.958 0.122
斜坡（可以）
极性 tEC.trend 东南 df lower.CL upper.CL
否定 87205.18 20992.98 4 28919.33 145491.0
正面 145338.35 20992.98 4 87052.51 203624.2
对比 估计 东南 df t.ratio p.value
负 - 正 -58133.17 29688.55 4 -1.958101 0.12182
拦截（问题）
极性 lsmean 东南 df lower.CL upper.CL
否定 711.6652 22.2867 4 649.7874 773.543
正面 1117.0787 22.2867 4 1055.2009 1178.957
对比 估计 东南 df t.ratio p.value
负 - 正 -405.4135 31.51816 4 -12.86285 0.000211

emmeans 计算的截距与 lm 计算的截距不同。我认为问题在于模型没有为 EC = 0 定义。但我不确定。

解决方法

你所说的拦截不是；它们是 tEC 平均值处的模型预测。如果您想要拦截，请改用：

m.intercept <- lsmeans(model,"Polarity",at = list(tEC = 0))

您可以通过以下方式判断正在使用的参考水平

ref_grid(model)     # or   str(m.intercept)

请注意，这里拟合的模型由两条斜率不同的线组成；因此，预测之间的差异取决于 tEC 的值。因此，我强烈建议反对测试截距的比较；这些是 tEC 值的预测，如您所说，甚至无法发生。相反，尽量不要像数学家一样，对 tEC 的几个有代表性的值进行比较，例如，

LSMs <- lsmeans(model,at = list(tEC = c(0.001,0.003,0.005)))
pairs(LSMs,by = tEC)

您还可以轻松地将拟合线可视化：

emmip(model,Polarity ~ tEC,cov.reduce = range)
contrastemmeansinterceptrr