效果与预测功能

Question

我同时试图理解 R 的 predict() 函数和“效果”包 effect() 函数。本质上，我正在运行回归以测试两个二分 IV 在 DV 上的交互作用，同时控制两个连续协变量。在我的实际数据集中，交互很重要，所以现在我想绘制交互。因为我的模型中有协变量，所以我应该在控制这些其他变量（即 SPSS 中的估计边际均值）之后绘制均值。我之前没有在 R 中这样做过，在搜索时我开始期望我应该能够获得我需要的值，用于使用 effect() 或 predict() 函数进行绘图。因此，我尝试在随机生成的数据集上使用每个：

> set.seed(100)
> test <- data.frame(iv1 = factor(round(rnorm(200, mean=.5, sd=.25), 0), levels=c(0,1), labels=c("A","B")), iv2 = factor(round(rnorm(200, mean=.5, sd=.25), 0), levels=c(0,1), labels=c("C","D")), cv1 = rnorm(200, mean=4, sd=1), cv2 = rnorm(200, mean=3, sd=1), dv = rnorm(200, mean=5, sd=1))
> mod <- lm(dv ~ cv1 + cv2 + iv1*iv2, data = test)
> new <- with(test, expand.grid(iv1 = levels(iv1), iv2 = levels(iv2), cv1 = mean(cv1), cv2 = mean(cv2)))
> test$pv <- predict(mod, newdata = new)

> tapply(test$pv, list(test$iv1, test$iv2), mean)
         C        D
A 5.076842 5.086218
B 5.025614 5.065399

> effect("iv1:iv2", mod)

 iv1*iv2 effect
   iv2
iv1        C        D
  A 5.019391 5.167275
  B 5.216955 4.855195

因为我得到不同的结果，所以我将数据导出到 SPSS 并运行 ANOVA 做同样的事情并查看估计的边际均值 (EMMEANS)。这些与 R 中 effect() 给出的结果相同。

SPSS 语法：

DATASET ACTIVATE DataSet1.
RECODE iv1 iv2 ('A'=-1) ('B'=1) ('C'=-1) ('D'=1) INTO iv1_recode iv2_recode.
EXECUTE.

UNIANOVA dv BY iv1_recode iv2_recode WITH cv1 cv2
  /METHOD=SSTYPE(3)
  /INTERCEPT=INCLUDE
  /EMMEANS=TABLES(OVERALL) WITH(cv1=MEAN cv2=MEAN) 
  /EMMEANS=TABLES(iv1_recode) WITH(cv1=MEAN cv2=MEAN) 
  /EMMEANS=TABLES(iv2_recode) WITH(cv1=MEAN cv2=MEAN) 
  /EMMEANS=TABLES(iv1_recode*iv2_recode) WITH(cv1=MEAN cv2=MEAN) 
  /PRINT=DESCRIPTIVE
  /CRITERIA=ALPHA(.05)
  /DESIGN=cv1 cv2 iv1_recode iv2_recode iv1_recode*iv2_recode.

作为检查，EMMEANS 的 SPSS 输出显示：“模型中出现的协变量按以下值进行评估：cv1 = 3.996208827095569, cv2 = 3.052881951477868。”这些与我用于预测的协变量的值相同：

> new
  iv1 iv2      cv1      cv2
1   A   C 3.996209 3.052882
2   B   C 3.996209 3.052882
3   A   D 3.996209 3.052882
4   B   D 3.996209 3.052882

那么我有什么不明白的呢？还是我在这里做一些愚蠢的事情（一种明显的可能性）？这可能是我没有掌握估计的边际均值。

非常感谢任何帮助！

Answer 1

因此，对于模型本身获取结果以及将模型应用于观察到的数据，这里似乎有些混乱。这里出现了一个大问题

test$pv <- predict(mod, newdata = new)

这里，new 有 4 行，所以 predict(mod, newdata = new) 有值。运行就是这样

predict(mod, newdata = new)
#        1        2        3        4 
# 5.019391 5.216955 5.167275 4.855195 

并注意这些值如何匹配来自effect() 的结果。

当您将它们分配给 test$pv 时，该长度为 4 的向量会被回收，因此它最终会沿着 test data.frame 重复 50 次。而test 确实包含您观察到的数据，因此混合模型的理论预测和观察到的数据并不是一个真正的好主意。如果您真的想要每个观察的“真实”预测值，那么test$pv<-predict(mod) 将是正确的选择。但是，将总和超过test，这又是观察值`，

tapply(test$pv, list(test$iv1, test$iv2), mean)

将使用实际观察到的 cv1 和 cv2 的值，而不仅仅是协变量的总体平均值。

我们已经看到 effect() 使用协变量的平均值，但您也可以使用显式设置值

effect("iv1:iv2", mod, given.values=c(cv1=3.996209, cv2=3.052882))

如果你喜欢。